Polynomials  are  one  of  principal  tools  of  classical  numerical  analysis.  When  a  function 
needs  to  be  interpolated,  integrated,  differentiated,  etc.,  it  is  assumed  to  be  approximated 
by  a  polynomial  of  a  certain  fixed  order  (though  the  polynomial  is  almost  never  constructed 
explicitly),  and  a  treatment  appropriate  to  such  a  polynomial  is  applied.  We  introduce  anal¬ 
ogous  techniques  based  on  the  assumption  that  the  function  to  be  dealt  with  is  band-limited, 
and  use  the  well-developed  apparatus  of  Prolate  Spheroidal  Wave  Functions  to  construct 
quadratures,  interpolation  and  differentiation  formulae,  etc.  for  band-limited  functions. 
Since  band-limited  functions  are  often  encountered  in  physics,  engineering,  statistics,  etc. 
the  apparatus  we  introduce  appears  to  be  natural  in  many  environments.  Our  results  are 
illustrated  with  several  numerical  examples. 
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Prolate  Spheroidal  Wave  Functions,  Quadrature,  and 

Interpolation 


1  Introduction 

Numerical  quadrature  and  interpolation  are  a  well-developed  part  of  numerical  analysis; 
polynomials  are  the  classical  tool  for  the  design  of  such  schemes.  Conceptually  speaking, 
one  assumes  that  the  function  is  well-approximated  by  expressions  of  the  form 

(1) 

j=0 

with  reasonably  small  n,  and  designs  algorithms  that  are  effective  for  functions  of  the 
form  (1)  (needless  to  say,  one  almost  never  actually  computes  the  coefficients  {aj};  one 
only  uses  the  fact  of  their  existence).  Obviously,  the  polynomial  approach  is  only  effective 
for  functions  that  are  well- approximated  by  polynomials. 

When  one  has  to  handle  functions  that  are  well-behaved  on  the  whole  line  (for  ex¬ 
ample,  in  signal  processing),  polynomials  are  not  an  appropriate  tool.  In  such  cases, 
trigonometric  polynomials  are  used;  existing  tools  are  very  satisfactory  for  dealing  with 
functions  defined  and  well-behaved  on  the  whole  of  R,^.  Such  tools,  in  effect,  make  the 
assumption  that  the  functions  are  band-limited  or  nearly  so;  a  function  /  ;  R,  -4-  R,  is 
said  to  be  band-limited  if  there  exist  a  positive  real  c  and  a  function  cr  €  L^[— 1, 1]  such 
that 

f{x)  =  £^e'“*a(t)  dt.  (2) 

However,  in  many  cases,  we  are  confronted  with  band-limited  functions  defined  on  inter¬ 
vals  (or,  more  generally,  on  compact  regions  in  R").  Wave  phenomena  are  a  rich  source 
of  such  functions,  both  in  the  engineering  and  computational  contexts;  they  are  also 
encountered  in  fluid  dynamics,  signal  processing,  and  many  other  areas.  Often,  such 
functions  can  be  effectively  approximated  by  polynomials  via  standard  tools  of  classical 
analysis.  However,  even  when  such  approximations  are  feasible,  they  are  usually  not 
optimal.  Smooth  periodic  functions  are  a  good  illustration  of  this  observation:  while 
they  can  be  approximated  by  polynomials  (for  example,  via  Chebyshev  or  Legendre 
expansions),  they  are  more  efficiently  approximated  by  Fourier  expansions,  both  for  an¬ 
alytical  and  numerical  purposes.  It  would  appear  that  an  approach  explicitly  based  on 
trigonometric  polynomials  could  be  more  efficient  in  dealing  with  band-limited  functions. 

In  the  engineering  context,  such  an  apparatus  was  constructed  more  than  30  years 
ago  (see  [20]-[21],  [7]-[9]).  The  natural  tool  for  analyzing  band-limited  functions  on  R^  is 
the  Fourier  Transform,  unless  the  functions  are  periodic,  in  which  case  the  natural  tool  is 
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the  Fourier  Series.  The  authors  of  [20]- [21]  observe  that  for  the  analysis  of  band-limited 
functions  on  the  interval,  Prolate  Spheroidal  Wave  Functions  are  likewise  a  natural  ap¬ 
proach.  The  authors  also  construct  a  multidimensional  version  of  the  theory,  though 
their  apparatus  is  only  complete  for  the  case  of  spherical  regions. 

The  present  paper  constructs  tools  for  the  use  of  the  approach  of  [20]-[21]  in  the 
modern  computational  environment.  We  construct  a  class  of  quadratures  for  band- 
limited  functions  that  closely  parallel  the  Gaussian  quadratures  for  polynomials.  The 
nodes  are  very  close  to  being  roots  of  appropriately  chosen  Prolate  Spheroidal  Wave 
Functions,  the  resulting  quadratures  are  stable,  and  all  weights  are  positive.  As  in  the 
case  of  polynomials,  there  are  interpolation,  differentiation  and  indefinite  integration 
schemes  associated  with  the  obtained  quadratures,  exact  on  certain  classes  of  band- 
limited  functions.  These  procedures  are  the  main  tools  necessary  for  the  numerical  use 
of  spectral  discretizations  based  on  Prolate  Spheroidal  Wave  Functions,  instead  of  on 
the  usual  polynomial  bases.  When  dealing  with  band-limited  functions,  the  number  of 
nodes  required  by  these  procedures  to  obtain  a  prescribed  accuracy  is  much  less  than 
that  required  by  their  pol>'nomial-based  counterparts.  An  additional  bonus  is  the  fact 
that  the  condition  number  of  differentiation  of  prolate  spheroidal  wave  functions  is  less 
than  that  of  differentiation  of  the  usual  polynomial  basis  functions  (see  Section  8  below). 

This  paper  is  organized  as  follows.  Section  2  summarizes  various  standard  mathemat¬ 
ical  facts  used  in  the  remainder  of  the  paper.  Section  3  contains  derivations  of  various 
results  used  in  the  algorithms  described  in  later  sections.  Section  4  describes  algorithms 
for  evaluation  of  prolate  spheroidal  wave  functions  and  associated  eigenvalues.  Section  5 
describes  a  construction  of  quadratures  for  band-limited  functions.  Section  6  describes 
an  alternative  approach  to  arriving  at  such  quadratures;  it  shows  that  roots  of  appropri¬ 
ately  chosen  prolate  spheroidal  wave  functions  can  serve  as  quadrature  nodes.  Section  7 
analyzes  the  use  of  prolate  spheroidal  wave  functions  for  interpolation.  Section  8  con¬ 
tains  results  of  our  numerical  experiments  with  quadratures  and  interpolation.  Section  9 
contains  a  number  of  miscellaneous  properties  of  prolate  spheroidal  wave  functions,  and 
Section  10  contains  generalizations  and  conclusions. 

2  Mathematical  Preliminaries 

As  a  matter  of  convention,  in  this  paper  the  norm  of  a  function  is,  unless  stated  otherwise, 
its  norm: 

ll/ll  =  j  dx.  (3) 
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2.1  Chebyshev  systems 


Definition  2.1  A  sequence  of  functions  (f)i, . . . ,  (pn  will  be  referred  to  as  a  Chebyshev 
system  on  the  interval  [a,  b]  if  each  of  them  is  continuous  and  the  determinant 


4>i{xi)  •••  P\{Xn) 

^71(2^1)  ■  ■  ■  Pn{Xji) 

is  nonzero  for  any  sequence  of  points  Xi,. .  .,Xn  such  that  a  <  Xi  <  X2  ■  ■  ■  <  Xn  <  b. 


(4) 


An  alternate  definition  of  a  Chebyshev  system  is  that  any  linear  combination  of  the 
functions  with  nonzero  coefficients  must  have  fewer  than  n  zeros. 

Examples  of  Chebyshev  and  extended  Chebyshev  systems  include  the  following  (ad¬ 
ditional  examples  can  be  found  in  [11]). 

Example  2.1  The  powers  l,a:, form  an  extended  Chebyshev  system  on  the 
interval  {—00,00). 

Example  2.2  The  exponentials  form  an  extended  Chebyshev  sys¬ 

tem  for  any  Ai, . . . ,  An  >  0  on  the  interval  [0, 00). 

Example  2.3  The  functions  1,  cos  x,  sin  x,  cos  2x,  sin  2a;, ... ,  cos  nx,  sin  nx  form  a  Cheby¬ 
shev  system  on  the  interval  [0, 27r]. 


2.2  Generalized  Gaussian  quadratures 

A  quadrature  rule  is  an  expression  of  the  form 

n 

YlwjPixj),  (5) 

j=i 

where  the  points  Xj  G  E  and  coefficients  Wj  G  E  are  referred  to  as  the  nodes  and  weights 
of  the  quadrature,  respectively.  They  serve  as  approximations  to  integrals  of  the  form 

rb 

/  <f){x)u{x)  dx,  (6) 

Ja 

with  u!  being  an  integrable  non-negative  function. 

Quadratures  are  typically  chosen  so  that  the  quadrature  (5)  is  equal  to  the  desired 
integral  (6)  for  some  set  of  functions,  commonly  polynomials  of  some  fixed  order.  Of 
these,  the  classical  Gaussian  quadrature  rules  consist  of  n  nodes  and  integrate  polynomi¬ 
als  of  order  2n  —  1  exactly.  In  [13],  the  notion  of  a  Gaussian  quadrature  was  generalized 
as  follows: 
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Definition  2.2  A  quadrature  formula  will  be  referred  to  as  Gaussian  with  respect  to  a 
set  of  2n  functions  (j)i,. . . ,  ^2n  :  [o,  — >■  Et  and  a  weight  function  u  :  [a,  6]  R"*",  if  it 

consists  of  n  weights  and  nodes,  and  integrates  the  functions  (l)i  exactly  with  the  weight 
function  u  for  all  i  =  1, . . . ,  2n.  The  weights  and  nodes  of  a  Gaussian  quadrature  will  he 
referred  to  as  Gaussian  weights  and  nodes  respectively. 

The  following  theorem  appears  to  be  due  to  Markov  [14,  15];  proofs  of  it  can  also  be 
found  in  [12]  and  [11]  (in  a  somewhat  different  form). 

Theorem  2.1  Suppose  that  the  functions  4)1, . . . ,  f)2n  ■  [o,  &]  R  form  a  Chebyshev 
system  on  [a,  6].  Suppose  in  addition  that  a;  :  [a,  6]  -4  R  is  a  non-negative  integrable 
function  [a,  b]  -4  R.  Then  there  exists  a  unique  Gaussian  quadrature  for  the  functions 
4>i}  -  ■  ■■>  4>2n  on  [a,  6]  with  respect  to  the  weight  function  co.  The  weights  of  this  quadrature 
are  positive. 

While  the  existence  of  Generalized  Gaussian  Quadratures  was  observed  more  than 
100  years  ago,  the  constructions  found  in  [14,  15],  [6,  12],  [10,  11]  do  not  easily  yield 
numerical  algorithms  for  the  design  of  such  quadrature  formulae;  such  algorithms  have 
been  constructed  recently  (see  [13,  25,  2]). 

Remark  2.1  It  might  be  worthwhile  to  observe  here  that  when  a  Generalized  Gaussian 
quadrature  is  to  be  constructed,  the  determination  of  its  nodes  tends  to  be  the  critical 
step  (though  the  procedure  of  [13,  25,  2]  determines  the  nodes  and  weights  simultane¬ 
ously).  Indeed,  once  the  nodes  xi,X2,  ■ . .  ,Xn  have  been  found,  the  weights  wi,W2, . . .  ,Wn 
can  be  determined  easily  as  the  solution  of  the  n  x  n  system  of  linear  equations 

■  4>i{xj)  =  /  4)i{x)  dx,  (7) 

■'O' 

with  i  =  1,2, . . .  ,n. 

2.3  Legendre  Polynomials 

In  agreement  with  standard  practice,  we  will  be  denoting  by  the  classical  Legendre 
polynomials,  defined  by  the  three-term  recursion 

Pn+i  (a:)  =  X  Pn  (a;)  -  (a:) ,  (8) 

with  the  initial  conditions 

Po{x)  =  1,  (9) 

Pi  (re)  =  x; 


4 


as  is  well-known 


n(l)  =  1  (10) 

for  all  A;  =  0, 1, 2, . . and  each  of  the  polynomials  satisfies  the  differential  equation 

(11) 

The  polynomials  defined  by  the  formulae  (8), (9)  are  orthogonal  on  the  interval  [-1,1]; 
however,  they  are  not  orthonormal,  since  for  each  n  >  0, 

dx  = (12) 

the  normalized  version  of  the  Legendre  polynomials  will  be  denoted  by  7^,  so  that 

Pn{x)  =  Pn{x)  •  ^n  +  1/2.  (13) 

The  following  lemma  follows  immediately  from  the  Cauchy-Schwartz  inequality  and  from 
the  orthogonality  of  the  Legendre  polynomials  on  the  interval  [—1,1]: 


Lemma  2.2  For  all  integer  k  >  n, 

r-l  .  _  I  /  2 


j  Pn{x)  dx 
For  all  integer  0  <  k  <n, 


< 


k+1 


X^'Pnix) 


dx 


=  0. 


(14) 


(15) 


2.4  Convolutional  Volterra  Equations 

A  convolutional  Volterra  equation  of  the  second  kind  is  an  expression  of  the  form 

rx 

V?(x)  =  /  K (x  -  t)  (p{t)  dt  -t-  a{x)  (16) 

Ja 

where  a,  b  are  a  pair  of  numbers  such  that  a  <  b,  the  functions  a,  K  :  [a,  b]  ^  are 
square-integrable,  and  y?  :  [a,  6]  ->  C  is  the  function  to  be  determined.  Proofs  of  the 
following  theorem  can  be  found  in  [4],  as  well  as  in  many  other  sources. 

Theorem  2.3  The  equation  (16)  always  has  a  unique  solution  on  the  interval  [a,  5].  If 
both  functions  K,  a  are  k  times  continuously  differentiable,  the  solution  (p  is  also  k  times 
continuously  differentiable. 
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2.5  Prolate  Spheroidal  Wave  Functions 

In  this  subsection,  we  summarize  certain  facts  about  the  Prolate  Spheroidal  Wave  Func¬ 
tions.  Unless  stated  otherwise,  all  these  facts  can  be  found  in  [20,  17], 

Given  a  real  c>  0,  we  will  denote  by  Fc  the  operator  L^[-l,  1]  ^  L^-l,  1]  defined 
by  the  formula 

F,{^){x)  =  (17) 

Obviously,  Fc  is  compact;  we  will  denote  by  Ao,  Ai, . . . ,  A„, . . .  the  eigenvalues  of  Fc 
ordered  so  that  |Aj_i|  >  |A_j|  for  all  natural  j.  For  each  non-negative  integer  j,  we  will 
denote  by  ipj  the  eigenfunctions  corresponding  to  Xj,  so  that 

Xjipjix)  =  j  dt,  (18) 

for  all  X  G  [—1, 1];  we  adopt  the  convention  that  the  functions  are  normalized  such  that 
ll'0illL2[_i,i]  =  1,  for  all  j}  The  following  theorem  is  a  combination  of  several  lemmas 
from  [20], [6], [11]. 

Theorem  2.4  For  any  positive  real  c,  the  eigenfunctions  '0o,'0i)  •  •  •  >  of  the  operator  Fc 
are  purely  real,  are  orthonormal,  and  are  complete  in  L^[— 1,1].  The  even-numbered 
eigenfunctions  are  even,  and  the  odd-numbered  ones  are  odd.  All  eigenvalues  of  Fc 
are  non-zero  and  simple;  the  even-numbered  eigenvalues  are  purely  real,  and  the  odd- 
numbered  ones  are  purely  imaginary;  in  particular,  Xj  =  P\Xj\.  The  functions  il^i  consti¬ 
tute  a  Chebychev  system  on  the  interval  [—1, 1];  in  particular,  the  function  ipi  has  exactly 
i  zeroes  on  that  interval,  for  any  i  =  0, 1, . . . ,. 

We  will  define  the  self-adjoint  operator  Qc :  T^[-l,  1]  T^[-l,  1]  by  the  formula 

=  i  (19) 

TT  X  —  t 

a  simple  calculation  shows  that 

0.  =  ^  •  J';  ■  Fc,  (20) 

that  Qc  has  the  same  eigenfunctions  as  Fc,  and  that  the  j-th  (in  descending  order) 
eigenvalue  p,j  of  Qc  is  connected  with  Xj  by  the  formula 

ho  ~  (21) 

^This  convention  differs  from  that  used  in  [20];  however,  the  present  paper  is  concerned  almost 
exclusively  with  approximation  of  functions  on  [—1, 1],  and  in  that  context,  the  convention  that  the 
functions  have  unit  norm  on  that  interval  is  by  far  the  most  convenient. 
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The  operator  Qc  is  obviously  closely  related  to  the  operator  Pc :  L^[-oo,  oo]  [-00, 00] 
defined  by  the  formula 


PM  =  -■  f 

TT  J- 


1  f°°  sin{c  •  {x  —  t)) 


X  —  t 


(22) 


which,  as  is  well  known,  is  the  orthogonal  projection  operator  onto  the  space  of  functions 
of  band  limit  c  on  (—00, 00). 

For  large  c,  the  spectrum  of  Qc  consists  of  three  parts:  about  2c/7r  eigenvalues  that 
are  very  close  to  1,  followed  by  order  log(c)  eigenvalues  which  decay  exponentially  from  1 
to  nearly  0;  the  remaining  eigenvalues  are  all  very  close  to  zero.  The  following  theorem, 
proven  (in  a  slightly  different  form)  in  [19],  describes  the  spectrum  of  Qc  more  precisely. 


Theorem  2.5  For  any  positive  real  c  and  0  <  a  <  1  the  number  N  of  eigenvalues  of 
the  operator  Qc  that  are  greater  than  a  satisfies  the  inequality 

7  +  (^  log  ^7^)  log(c)  -  10  •  log(c)  <N<  (23) 

7  +  (7  los(c)  +  10  •  log(c)- 

By  a  remarkable  coincidence,  the  eigenfunctions  -  •  ,'ipnOf  the  operator  Qc  turn 
out  to  be  the  Prolate  Spheroidal  Wave  functions,  well-known  from  classical  Mathematical 
Physics  (see,  for  example,  [16]).  The  following  theorem  formalizes  this  statement;  it  is 
proven  in  a  considerably  more  general  form  in  [21]. 


Theorem  2.6  For  any  c  >  0,  there  exists  a  strictly  increasing  sequence  of  positive  real 
numbers  Xo,Xi,--  ■  such  that  for  each  j  >  0,  the  differential  equation 

(1  -  7)  't]j"{x)  -  2x  ip'{x)  -I-  (xj  -  x^)  'ip(x)  =  0  (24) 

has  a  solution  that  is  continuous  on  the  interval  [—1, 1].  For  each  j  >  0,  the  function  'tpj 
(defined  in  Theorem  2.4)  is  the  solution  of  (24). 


3  Analytical  Apparatus 

3.1  Prolate  Series 

Since  the  functions  ^0,  ^1,  •  •  • ,  •  •  •  are  a  complete  orthonormal  basis  in  L^[— 1, 1],  any 

formula  for  the  inner  product  of  prolate  spheroidal  wave  functions  with  another  function 
/  is  also  a  formula  for  the  coeflBcients  of  an  expansion  of  /  into  prolate  spheroidal  func¬ 
tions  (which  we  will  refer  to  as  the  prolate  expansion  of  /).  Thus  the  following  theorem 
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provides  the  coefficients  of  the  prolate  expansion  of  the  derivative  of  a  prolate  spheroidal 
function,  and  also  the  coefficients  of  the  prolate  expansion  of  a  prolate  spheroidal  wave 
function  multiplied  by  x.  Those  coefficients  are  also  the  entries  of  the  matrix  for  differen¬ 
tiation  of  a  prolate  expansion  (producing  another  prolate  expansion),  and  the  entries  of 
the  matrix  for  multiplication  of  a  prolate  expansion  by  x,  respectively.  (These  formulae 
are  not,  however,  suitable  for  producing  such  matrices  numerically,  since  in  many  cases 
they  exhibit  catastrophic  cancellation.) 


Theorem  3.1  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  If  m  =  n  (mod  2),  then 


1  ^'<Pn{x)'>Pmix)  dx  = 

J  ^  X  '(pm{x)  dx  =  0. 

(25) 

Ifm^n  (mod 2),  then 

j  ^tp'^ix'j^rnix)  dx 

2  \2 

(26) 

J  ^X'tpn{o^)i^m{x)  dx 

(27) 

Proof.  Since  the  functions  -ipj  are  alternately  even  and  odd,  (25)  is  obvious.  In  order  to 
prove  (26),  we  start  with  the  identity 

K'ipn  -  j  ^  ipn{t)  dt  (28) 

(see  (18)  in  Subsection  2.5).  Differentiating  (28)  with  respect  to  x,  we  obtain 

Xn'ip'M  =  ic  f  V’n  (i)  dt.  (29) 

Projecting  both  sides  of  (29)  on  ifrn  and  using  the  identity  (28)  (with  n  replaced  with 
m)  again,  we  have 

K  Ipmix)  dx 

=  ic  J  1pm{x)  J  te^‘^^'lfn{t)  dtdx 
=  ic  j  tipnit)  J  e*“‘'0m(2;)  dxdt 

=  icXm  j  ^t'llJn{t)'lpm{t)  dt.  (30) 
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Obviously,  the  above  calculation  can  be  repeated  with  m  and  n  exchanged,  yielding  the 
identity 

Xm  dx  =  i  C  Xn  J  1pn{t)  dt\  (31) 

combining  (30)  with  (31),  we  have 

/  ^n{x)  dx  =  ^  f  xl}m{x)'tp'n{x)  dx.  (32) 

^-1  J-l 

On  the  other  hand,  integrating  the  left  side  of  (32)  by  parts,  we  have 
j^^^!^{x) 'tpnix)  dx 

=  tprnil)  ^n{l)  “  '0m(-l)  ^n(-l)  “  J_^'(p'r,ix)  dx.  (33) 

Since  m^n  (mod 2),  we  rewrite  (33)  as 

'fpnix)  dx 

=  2  ^to(I)  ^„(1)  -  Tp'^{x)  'Ipm.{x)  dx.  (34) 

Now,  combining  (32)  and  (34)  and  rearranging  terms,  we  get 

/  ^n(^)  i’mix)  dx  =  ^^(1)  t/;„(l)  .  (35) 

•^-1  +  Xi 

Substituting  (30)  into  (35),  we  get 


J  ^X'll>n{x)lprn{x)  dx 

=  /  i^nix)  i^rnix)  dx 

*'“1 


ic  Xm  J-l 

1  A„  2X1 
ic  XjTfi 

2  XjnXji 

ic  A2  +  A2 

771  '  71 


^m(l)'0n(l) 


(36) 

□ 


The  following  corollary,  which  is  an  immediate  consequence  of  (32),  finds  use  in  the 
numerical  evaluation  of  the  eigenvalues  {Aj}: 
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Corollary  3.2  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  Ifm^n  (mod 2),  then 


A2  f  tpmix)  dx 

''m  _  J-l _ 


(37) 


3.2  Decay  of  Legendre  Coefficients  of  Prolate  Spheroidal  Wave- 
functions 

Since  each  of  the  functions  'ihj  is  analytic  on  C,  on  the  interval  [-1, 1]  it  can  be  expanded 
in  a  Legendre  series  of  the  form 

00 

=  Y.  PkPk{x),  (38) 

k=0 

with  the  coefficients  Pk  decaying  superalgebraically;  the  following  two  theorems  establish 
bounds  for  the  decay  rate. 


Lemma  3.3  Let  P„(x)  be  the  n-th  normalized  Legendre  polynomial  (defined  in  (13)). 
Then  for  any  real  a, 

j  e^°'^Pn{x)dx 

oo  fl  _  00  .1 

=  Y  x^'^  Pn{x)  dx i  Y  l^k  Pn{x)  dx.  (39) 

k=ko  k=ko 


where 


Cti  =  (-1)' 

A  =  (-1)* 


(2*:)!’ 


(2A:  +  1)!’ 


ko  =  [n/2j . 

Furthermore,  for  all  integer  m  >  [e  •  |a|J  +  1, 


/•l  _  rn-l  .1 

/  e^^=^Pn{x)dx-  Y.  ak 

77^1  1 

-iT.h 

/i.io 


X  Pn{x)  dx 


x'^  Pn{x)  dx 


(40) 

(41) 

(42) 
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(43) 


In  particular,  if 


n>2  (Le-  |a|J  +1)  , 


(44) 


then 


e‘“  P,{x) 


dx 


jN  n-l 

<'2, 


(45) 


Proof.  The  formula  (39)  follows  immediately  from  Lemma  2.2  and  Taylor’s  expansion 
of  In  order  to  prove  (43),  we  assume  that  m  is  an  integer  such  that 


m>  [e-  |o|J  +  1 . 

Introducing  the  notation 

Pn{x)  dx-\-i'Y^  13k  /  Pn{x)  dx, 

k=m  k=m 

we  immediately  observe  that,  due  to  Lemma  2.2  and  the  triangle  inequality. 


l-R™!  <  E 

k=2m. 


m 


<  E 


a 


kl  k  +  1 

\k 


k=2m 


k\  ■ 


Since  (46)  implies  that 

lal 


1 


2m  +  k  ^  2m  ^  2e  2  ’ 
for  all  integer  m,k  >  0,  we  rewrite  (48)  as 

12771 


1^41  < 


a 


(2m) ! 

12771 


7  1  1 

1+2+J  +  --. 


<  2 


a 


(46) 


(47) 


(48) 


(49) 


(2m) ! 


(50) 


and  obtain  (43)  immediately  using  Stirling’s  formula.  Finally,  we  obtain  (45)  by  choosing 
m  =  [e  •  |a|J  +  1 .  (51) 

□ 
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Theorem  3.4  Let  'ipm{x)  he  the  m-th  prolate  spheroidal  function  with  band  limit  c,  let 
Pk{x)  be  the  k-th  normalized  Legendre  polynomial  (defined  in  (18)),  and  let  Xm  be  the 
eigenvalue  which  corresponds  to  -ipruix)  (as  in  Theorem  2.4)-  Then  for  all  integer  m  >  0 
and  all  real  positive  c,  if 


A:  >  2  ([e-cj  +1)  , 

(52) 

then 

fi  _  1 

J  ^  i;rn{x)  Pk{x)  ^  ~  '  [2) 

(63) 

Moreover,  given  any  e  >  0,  if 

A:  >  2  (Le  •  cj  +  1)  +  log2  (j)  +  logs  (3^) , 

(54) 

then 

|y ^  ipm{x)Pi){x)  dx  <e. 

(55) 

Proof.  Obviously 

f  7pm{x)Pf:{x)  dx 

1/  —  1 

“  lAyy,!  \l-l  il-l  ^ 

<  \\  \  f  \'^m{x)\-  [  Pk{t)  dt  dx . 

\^m\  •'“1  J —  \ 

(56) 

Introducing  the  notation 


a  =  cx, 

and  remembering  that 

\'il}m{x)\dx  =  1, 

we  observe  that  the  combination  of  (56),  (57),  (58),  and  Lemma  3.3  implies  that 
j^^fiim{x)Pk{x)  dx\^ 

1  P 


Substituting  (54)  into  (53),  we  immediately  see  (55). 


(57) 

(58) 

(59) 
□ 
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4  Numerical  Evaluation  of  Prolate  Spheroidal  Wave- 
functions 

Both  the  classical  Bouwkamp  algorithm  (see,  for  example,  [1])  for  the  evaluation  of  the 
functions  ipj,  and  the  algorithm  presented  in  this  paper  for  the  same  task,  are  based  on 
the  expression  of  those  functions  as  a  Legendre  series  of  the  form 


AW  =  S  “WtW; 


(60) 


fc=0 


since  the  functions  'tpj  are  smooth,  the  coefficients  decay  superalgebraically  (with 
bounds  for  that  decay  being  given  in  Theorem  3.4).  Substituting  (60)  into  (24),  and 
using  (8)  and  (11),  we  obtain  the  well-known  three-term  recursion 


^ -f  2)(A:  4- 1)  2  ^  , 

•  C  •  Q!*+2  + 


{2k  +  3){2k  +  5) 


'  ,  2A:(A:4-1)-1  2  ^ 

k{k-l)  2 


(61) 


c  •  ak-2  =  0. 


{2k  -  3){2k  -  1) 

Combining  (61)  with  (13),  we  obtain  the  three-term  recursion 

{k  +  2){k  +  l) 


{2k  +  3)  ^  {2k +  5)  {2k +  1) 


■  PU2  + 


k{k  +  1)  + 


2k{k  -41)  -  1 


{2k  +  3){2k-l) 
k{k-l)  2 


-  Xj  •  + 


(62) 


{2k  -  l)^{2k-3){2k  +  l) 
for  the  coefficients  /3o,  /3i,. ..  of  the  expansion 

00 

=  J2^k-  'K{x)\ 


^  '  Pk-2  ~  0 


(63) 


k=0 


for  each  j  =  0, 1, 2, . . . ,  we  will  denote  by  the  vector  in  f  defined  by  the  formula 

p^  =  {pi,pipi,...).  (64) 

The  following  theorem  restates  the  recursion  (62)  in  a  slightly  different  form. 
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Theorem  4.1  The  coefficients  Xi  the  eigenvalues  and  the  vectors  /?*  are  the  corre¬ 
sponding  eigenvectors  of  the  operator  — >■  represented  by  the  symmetric  matrix  A 
given  by  the  formulae 


Ak,k 

Ak,k+2 

Ak+2,k 


^  (2fc  +  3)(2A:-l) 

{k  +  2)(A:  +  1) 

{2k  +  3)^(2A:  +  l)(2A:  +  5) 
{k  +  2)(A;  +  1) 

(2k  +  3)\/(2«:  +  l)(2*:  +  5) 


(65) 

(66) 

(67) 


for  all  A:  =  0, 1, 2, . . with  the  remainder  of  the  entries  of  the  matrix  being  zero. 


In  other  words,  the  recursion  (62)  can  be  rewritten  in  the  form 

(A-X,-/)(3^)=0,  (68) 

where  A  is  separable  into  two  symmetric  tridiagonal  matrices  Aeven  and  Aodd,  the  first 
consisting  of  the  elements  of  .4  with  even-numbered  rows  and  columns  and  the  second 
consisting  of  the  elements  of  .4  with  odd-numbered  rows  and  columns.  While  these  two 
matrices  are  infinite,  and  their  entries  do  not  decay  much  with  increasing  row  or  column 
number,  the  eigenvectors  {3^}  of  interest  (those  corresponding  to  the  first  m  prolate 
spheroidal  functions)  lie  almost  entirely  in  the  leading  rows  and  columns  of  the  matrices 
(as  shown  by  Theorem  3.4).  Thus  the  evaluation  of  prolate  spheroidal  functions  can  be 
performed  by  the  following  procedure: 


•  1.  Generate  the  leading  k  rows  and  columns  of  A,  where  k  is  given  by  (54). 

•  2.  Split  the  generated  portion  of  A  into  A^ven  and  A^dd,  and  use  a  solver  for  the 
symmetric  tridiagonal  eigenproblem  (such  as  that  in  LAPACK)  to  compute  their 
eigenvectors  {/?^}  and  eigenvalues  {xj}. 

•  3.  Use  the  obtained  values  of  the  coefficients  /?o,  /?2,  •  •  •  in  the  expansion  (63)  to 

evaluate  the  function  ifj  at  arbitrary  points  on  the  interval  [—1, 1]. 

Obviously  steps  1  and  2  can  be  performed  as  a  precomputation,  for  any  given  value  of 
c.  As  a  numerical  diagonalization  of  a  positive  definite  tridiagonal  matrix  with  well- 
separated  eigenvalues,  this  precomputation  stage  is  numerically  robust  and  efficient, 
requiring  0{cm)  operations  to  construct  the  Legendre  expansions  of  the  form  (64)  for  the 
first  m  prolate  spheroidal  functions;  each  subsequent  evaluation  of  a  prolate  spheroidal 
function  takes  0(c)  operations. 
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4.1  Numerical  Evaluation  of  Eigenvalues 

Although  the  above  algorithm  for  the  evaluation  of  prolate  spheroidal  wave  functions  also 
produces  the  eigenvalues  {xj}  of  the  differential  operator  (24),  it  does  not  produce  the 
eigenvalues  {Xj}  of  the  integral  operator  Fc  (defined  in  (17)).  Some  of  those  eigenvalues 
can  be  computed  using  the  formula 

Xj'ipjix)  =  f  dt,  (69) 

J  —  X 

evaluating  the  integral  on  the  right  hand  side  numerically;  however,  that  evaluation 
obviously  has  a  condition  number  of  about  l/Xj,  and  is  thus  inappropriate  for  computing 
small  Xj.  A  well-conditioned  procedure  is  as  follows: 

•  1.  Use  (69)  to  calculate  Aq,  evaluating  the  right  hand  side  numerically,  and  with 
a;  =  0  (so  that  ipo{x)  is  not  small). 

•  2.  Use  the  calculated  Ao,  together  with  Corollary  3.2,  to  compute  the  absolute  val¬ 
ues  |Aj|,  for  j  =  1, 2, ...,  m,  computing  each  |Aj|  from  \Xj-i\  (and  again,  evaluating 
the  required  integrals  numerically). 

•  3.  Use  the  fact  that  Xj  =  P\Xj\  (see  Theorem  2.4)  to  finish  the  computation. 

5  Quadratures  for  Band-Limited  Functions 

Since  the  prolate  spheroidal  wave  functions  tpo,ipi, . . .  ■  ■  ■  constitute  a  complete  or¬ 

thonormal  basis  in  L^[—l,  1]  (see  Theorem  2.4), 

=  f;  (  r  e-- ^.(r)  dr)  i^j(t),  (70) 

j=o  v-i  / 

for  all  z,  t  €  [—1, 1];  substituting  (18)  into  (70)  yields 

00 

j=0 

Thus  if  a  quadrature  integrates  exactly  the  first  n  eigenfunctions,  that  is,  if 


(71) 


for  all  j  =  0, 1, . . . ,  n  —  1,  then  the  error  of  the  quadrature  when  applied  to  a  function 
f{x)  =  with  a  e  [-1,1],  is  given  by 

m  I 

T  dx 


=  (a)  i’jixk)]  -  [  (  £  Xj  (a)  (rr)  ]  dx 

k=l  \j=o  J  “'-1  \j=0  ) 


I  \  I ,  ^  \ 

=  S  IT  -  /_  I]  Aj  'ipj{a)  'tpjix)  dx.  (73) 

A:=l  \j=n  /  ^  \j=n  / 

Due  to  the  orthonormality  of  the  functions  {V’j}, 


Y,Xj'tpj{a)ijj{x) 


3—n 


,,  E  l^iP- 

\j=n 


(74) 


From  (74),  it  is  obvious  that  the  error  of  integration  (73)  is  of  roughly  the  same  mag¬ 
nitude  as  An,  provided  that  n  is  in  the  range  where  the  eigenvalues  {Aj}  are  decreasing 
exponentially  (as  is  the  case  for  quadratures  of  any  useful  accuracy;  see  Theorem  2.5) 
and  provided  in  addition  that  the  weights  {wk}  are  not  large. 

Now,  the  existence  of  an  n/2-point  quadrature  that  is  exact  for  the  first  n  Prolate 
Spheriodal  Wave  functions  follows  from  the  combination  of  Theorems  2.1,  2.4;  an  al¬ 
gorithm  for  the  numerical  evaluation  of  nodes  and  weights  of  such  quadratures  can  be 
found  in  [2].  An  alternative  procedure  for  the  construction  of  quadrature  formulae  for 
band-limited  functions  (leading  to  slightly  different  nodes  and  weights)  is  described  in 
the  following  section;  a  numerical  comparison  of  the  two  can  be  found  in  Section  8  below. 


Remark  5.1  The  above  text  considers  only  the  error  of  integration  of  a  single  exponen¬ 
tial.  For  a  band-limited  function  g  :  [-1, 1]  -4  C  given  by  the  formula 

g{x)  =  (75) 

for  some  function  G  :  [— 1, 1]  -4-  C,  the  error  is  obviously  bounded  by  the  formula 

<£-||C?ll,  (76) 

where  e  is  the  maximum  error  of  integration  (73)  of  a  single  exponential,  for  any  t  € 
[— 1, 1].  While  |1G||  might  be  much  larger  than  ||^||[_i,i]  (as  it  is  if,  for  instance,  g  =  i’sQ.n), 
if  the  same  equation  (75)  is  used  to  extend  g  to  the  rest  of  the  real  line,  then  by  Parseval’s 
formula  ||G||  =  ||5||(-oo,oo);  that  is  to  say,  although  the  error  of  such  a  quadrature  when 
applied  to  a  band-limited  function  is  not  bounded  proportional  to  the  norm  of  that 
function  on  the  interval  of  integration,  it  is  bounded  proportional  to  the  norm  of  that 
function  on  the  entire  real  line. 


22  Wkg{xk)  -  /  g{x)  dx 

fc=i 
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6  Quadrature  Nodes  from  Roots  of  Prolate  Func¬ 
tions 

An  alternative  to  the  approach  of  the  previous  section  is  to  use  roots  of  appropriate 
prolate  spheroidal  wave  functions  as  quadrature  nodes,  with  the  weights  determined  via 
the  procedure  described  in  Remark  2.1.  The  following  theorems  provide  a  basis  for  this; 
numerically  (see  Section  8)  the  resulting  quadrature  nodes  tend  to  be  inferior  to  those 
produced  by  the  optimization  scheme  of  [13,  25,  2];  however,  they  are  useful  as  starting 
points  for  that  scheme,  or  as  somewhat  less  efficient  nodes  which  can  be  computed  much 
more  quickly. 

6.1  Euclid  Division  Algorithm  for  Band-Limited  Functions 

The  following  two  theorems  constitute  a  straightforward  extension  to  band-limited  func¬ 
tions  of  Euclid’s  division  algorithm  for  polynomials'.  Their  proofs  are  quite  simple,  and 
are  provided  here  for  completeness,  since  the  author  failed  to  find  them  in  the  literature. 

Theorem  6.1  Suppose  that  a,  :  [0, 1]  -4  C  are  a  pair  of  c^— functions  such  that 


<^(1)  ^  0,  (77) 

c  is  a  positive  real  number,  and  the  functions  f,  p  are  defined  by  the  formulae 

f{x)=l'  a{t)e‘^'‘'''*  dt,  (78) 

J  0 

p{x)  =  f  ip(t)  6*“’*  dt.  (79) 

J  0 

Then  there  exist  two  c^ -functions  77,^  :  [0, 1]  C  such  that 

f{x)  =  p{x)  q{x)  -h  r{x)  (80) 

for  all  3;  €  R,  with  the  functions  q,  r  :  [0, 1]  -4  R  defined  by  the  formulae 

q{x)  =  j  ?7(t)  dt,  (81) 

J  0 

r{x)  =  l\{t)e^^Ut.  (82) 
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Figure  1:  The  split  of  integration  range  that  yields  (85) 


Proof. 

Obviously,  for  any  functions  p,  g  given  by  (79),  (81), 

p{x)  q{x)  =  /  pit)  dt  ■  f  t)(t)  dr 

Jo  Jo 

=  drdt  (83) 

Defining  the  new  independent  variable  u  by  the  formula 

u  =  t  +  T,  (84) 

we  rewrite  (83)  as 

p{x)q{x)  =  f  [  ip{u—r)rj(T)  dr  du 

Jo  Jo 

+  [  f  <p{u—T)p{r)  dr  du  (85) 

Jl  Ju—l 

(see  Figure  1).  Substituting  (78),  (82),  and  (85)  into  (80),  we  get 

[  [  (p{u—T)  rjir)  drdu 

Jo  Jo 

+  (P{u-t)  77(t)  dr  du  +  ^{t)  dt 

Jl  Ju—l  Jo 

/■1/2  /•! 

=  /  a(t)  dt  +  /  a{t)  6^*“*  dt.  (86) 

^0  Jl/2 

Due  to  the  well  known  uniqueness  of  the  Fourier  Transform,  (86)  is  equivalent  to  two 
independent  equations: 

[  [  ip{u-r)r}{r)  dr  du-\-  [  ^it)e'‘^^  dt 

Jo  Jo  Jo 

/•1/2 

=  /  a{t)  dt,  (87) 

Jo 


18 


gicux  p  ^  ^(u-r)  77(t)  dr  du  =  ^  a{t)  dt. 


Now,  we  observe  that  (88)  does  not  contain  and  use  it  to  obtain  an  expression  for  r] 
as  a  function  of  cp,  a.  After  that,  we  will  view  (87)  as  an  expression  for  ^  via  (p,  a,  rj. 
From  (88)  and  the  uniqueness  of  the  Fourier  Transform,  we  obtain 


/  cp{u-T)r]{T)  dr  =  a{-), 
Ju—l  L 


for  all  u  €  [1, 2].  Introducing  the  new  variable  v  via  the  formula 

u  =  li  -  1,  (90) 

we  convert  (89)  into 

f  tp(^y+i-r)r}{T)  dr  =  (91) 

Jv  L 

which  is  a  Volterra  equation  of  the  first  kind  with  respect  to  77;  differentiating  (91)  with 
respect  to  u,  we  get 


-(p{l)r}{v)  +  (p’{v+l-T)r]{T)  dr  = 


which  is  a  Volterra  equation  of  the  second  kind.  Now,  the  existence  and  uniqueness 
of  the  solution  of  (92)  (and,  therefore,  of  (89)  and  (88))  follows  from  Theorem  2.3  of 
Section  2. 

With  77  defined  as  the  solution  of  (89),  we  use  (87)  together  with  the  uniqueness  of 
the  Fourier  Transform,  to  finally  obtain 


P>iu-T)  77(r)  dr, 


for  all  u  e  [0, 1]. 

□ 

The  following  theorem  is  a  consequence  of  the  preceding  one. 

Theorem  6.2  Suppose  that  a,  (p  :  [—1, 1]  — )•  C  are  a  pair  of  c^— functions  such  that 
#  0,  p>{l)  7^  0,  c  is  a  positive  real  number,  and  the  functions  f,p  are  defined  by 
the  formulae 


f{x)  =  j  ^  a{t)  dt, 
p{x)  =  J  ip{t)  dt. 
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Then  there  exist  two  -functions  rj,^  :  [—1, 1]  -4  C  such  that 

f{x)  =  p{x)  q(x)  +  r{x)  (96) 

for  all  a:  €  R,  with  the  functions  r  :  [— 1, 1]  -4  R  defined  by  the  formulae 

V{t)  dt,  (97) 

r{x)=  I  ^{t)e^^*dt.  (98) 

J  —  1 

Proof. 

Defining  the  functions  /^,  f-,p+,  p-,  by  the  formulae 

Mx)  =  f  a(t)  dt,  (99) 

Jo 

f.(x)  =  (100) 

P+(^)  =  [  dt,  (101) 

Jo 

P-{x)  =  f  >^(t)e'“*dt,  (102) 

J  —  1 

we  observe  that  for  all  x  €  R* , 

f{x)  =  U{x)  + f.{x),  (103) 

p(x)  =p+(x) +p_(x).  (104) 

Due  to  Theorem  6.1,  there  exist  such  77+,  77_,  that 

/+(x)  =  p+(x)  q^{x)  +  r+(x),  (105) 

/-  (^)  =  P-  (2;)  Q-  (x)  +  r_  (x) ,  (106) 

with  the  functions  g+,  g_,  r+,  r_  defined  by  the  formulae 

q+{x)=  f  77+(t)e*“*  dt,  (107) 

Jo 

q-{x)  =  j  ^  77_(t)  dt,  (108) 
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r+i^)=[  oft,  (109) 

Jo 

r_  (x)  =  1°  (t)  dt.  (110) 

Now,  defining  q,  by  the  formula 

q{x)  q-{x)  +  q+{x)  (111) 

for  all  x  G  [—1, 1],  we  have 

p{x)q{x)  =  {p.{x)  +  p+{x))  ■  {q-{x)  +  q+{x)) 

=  P+ix)q+{x)+P-{x)q-{x)+p-{x)q+{x) +p+{x)q-{x),  (112) 

and  we  define  r{x)  by  the  obvious  formula 

r{x)  =  r.{x)  +  r+{x)  -  {p.{x)q+{x)  +  p+{x)  q.{x)).  (113) 

□ 


6.2  Quadrature  nodes  from  the  division  theorem 

In  much  the  same  way  that  the  division  theorem  for  polynomials  can  be  used  to  provide 
a  constructive  proof  of  Gaussian  quadratures,  Theorem  6.2  provides  a  method  of  con¬ 
structing  generalized  Gaussian  quadratures  for  band-limited  functions.  The  method  is 
as  follows. 

To  construct  a  quadrature  for  functions  of  a  bandwidth  2c,  prolate  spheroidal  wave 
functions  corresponding  to  bandwidth  c  are  used.  (Thus  the  eigenvalues  {Aj}  and  eigen¬ 
functions  {'ipj}  are  in  this  section,  as  elsewhere  in  the  paper,  those  corresponding  to 
bandwidth  c).  The  following  theorem  provides  a  bound  of  the  error  of  a  quadrature 
whose  nodes  are  the  roots  of  the  n’th  prolate  function  when  applied  to  a  function 
/  which  satisfies  the  conditions  of  the  division  theorem,  in  terms  of  the  norms  of  the 
quotient  and  remainder  of  /  divided  by  'ipn- 

Theorem  6.3  Suppose  that  xi,X2,  ■ .  ■  ,Xn  G  K,  are  the  roots  of 'ipn-  Let  the  numbers 
wi,W2, . . . ,  tCn  €  R  6e  such  that 

n 

Y^Wk'tpjixk)  =  I  i^jix)  dx,  (114) 

k=l 
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for  all  j  =  0, 1, . . . ,  n  —  1.  Then  for  any  function  f  :  [—1, 1]  -4-  C  which  satisfies  the 
conditions  of  Theorem  6.2, 

'^kfixk)  -  [  fix)  dx 

A:=l 

00  /  m  \ 

<  |A»I  ■  Itell  +  Ilf II  ■  E  |A,I  ■  llv>iliL  •  2  +  E  .  (115) 

j=n  \  k=l  } 

where  the  functions  rj,^  :  [—1, 1]  -4  C  are  as  defined  in  Theorem  6.2. 

Proof.  Since  /  satisfies  the  conditions  of  Theorem  6.2,  there  exist  functions  q,r  : 
[— 1, 1]  -4  R,  defined  by  (97), (98)  such  that 

fix)  =  '(pnix)  qix)  +  r(a:).  (116) 

Then,  defining  the  error  of  integration  Ef  for  the  function  /  by 

Ef  =  '^kfixk)  -  f  fix)  dx 

k=l 


we  have 


n  -1 

=  Yl'fXkiifnixk)qixk) +rixk))  -  ii^nix)  qix)  +  rix))  dx 

k=l 

n  -1 

<  y^Wk  'tpnixk)  qixk)  -  /  fpnix)  qix)  dx 
k=l 

n  -1 

+  '^krixk)  -  /  rix)  dx 

l _ 1  v/— 1 


Since  the  nodes  {a;*}  are  the  roots  of 


y2'^k'ipnixk)qixk)  -  0. 

k=l 


Thus 


<  1/  ‘ipnix)qix)  dx  +  ^iOfcr(xfc)  -  [  rix)  dx 

J  tfnix)  qix)dx  =  J  fpnix)  j  r)it)  dt dx 

=  j  rjit)  J  'ipnix)e^‘^^  dxdt 


(118) 


(119) 


(120) 


(121) 
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Using  the  Cauchy-Schwartz  inequality  and  the  fact  that  the  function  'ipn  has  unit  norm, 
we  get  from  (121)  that 


Also, 


^„(x)  q{x)  dx 


<  |An|-|N|. 


n  .1 

Y,  Wkr{xk)  -  /  r{x) 

k=\ 


dx 


(122) 


(123) 


Substituting  (73)  into  (123),  and  using  the  Cauchy-Schwartz  inequality,  we  get 

n  -1 

Wkr{xk)  -  /  r{x)  dx 

k=l 

-1  f  m  /  oo  ' 

^  \k=l  \j=n  j 

-  j_(Y  Mi)  dx  I  dt 


\j=n 


<  ll^ll  •  E  M  ■  WMlo  •  2  +  X]  ||u;fcl|  . 

j=n  \  k=l  J 


Combining  (120),  (122),  and  (124),  we  get 


Ef  <  |A„i-M  +  Ilfll'El^li-WlIlL-  2  +  Ell“'*ll  • 

j=n  V  fc=l  / 


(124) 


(125) 

□ 


Remark  6.1  The  use  of  Theorem  6.3  for  the  construction  of  quadrature  rules  for  band- 
limited  functions  depends  on  the  fact  that  the  norms  of  the  band-limited  functions  q 
and  r  in  (116)  are  not  large,  compared  to  the  norm  of  /  (both  sets  of  norms  being  on 
[—00,00]).  Such  estimates  have  been  obtained  for  all  n  >  2c/7r  -I-  lOlog(c).  The  proofs 
are  quite  involved,  and  will  be  reported  at  a  later  date.  In  this  paper,  we  demonstrate 
the  performance  of  the  obtained  quadrature  formulae  numerically  (see  Section  8  below). 
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Remark  6.2  It  is  natural  to  view  (116)  as  an  analogue  for  band-limited  functions  of 
the  Euclid  division  theorem  for  polynomials.  However,  there  are  certain  differences.  In 
particular,  Theorem  6.1  admits  extensions  to  band-limited  functions  of  several  variables, 
while  the  classical  Euclid  algorithm  does  not.  Such  extensions  (together  with  several 
applications)  will  be  reported  at  a  later  date. 


7  Interpolation  via  Prolate  Spheroidal  Wavefunctions 

Interpolation  is  usually  performed  by  the  following  general  procedure;  assuming  that  the 
function  /  :  [a,  6]  C  to  be  interpolated  is  given  by  the  formula 

f{x)  =  Ci<f)i{x)  +  C2^2{x)  +  ...  +  Cn^nix),  (126) 

where  (t>i,(t>2,  ■  ■  ■  ,(i>n  '■  [a,  &]  C  are  a  fixed  sequence  of  functions  (often  polynomials), 
solve  an  n  X  n  linear  system  to  determine  the  coefficients  Ci,  C2, . . . ,  c„  from  the  values  of 
/  at  the  n  interpolation  nodes,  then  use  (126)  to  evaluate  /  wherever  needed.  As  is  well 
known,  if  /  is  well-approximated  by  a  linear  combination  of  the  interpolation  functions, 
and  if  the  linear  system  to  be  solved  is  well-conditioned,  then  this  procedure  is  accurate. 

As  shown  in  Section  5  in  the  context  of  quadratures,  a  linear  combination  of  the  first 
n  prolate  spheroidal  functions  for  a  band  limit  c  can  provide  a  good 

approximation  to  functions  of  the  form  with  t  G  [-1, 1]  (see  (71,74));  in  the  regime 
where  the  accuracy  is  numerically  useful,  the  error  is  of  the  same  order  of  magnitude  as 
|A„|.  This,  in  turn,  shows  that  they  provide  a  good  approximation  (in  the  same  sense  as 
in  Remark  5.1)  to  any  band-limited  function  of  band  limit  c.  Thus,  if  ^/’o,  V’l, . . . ,  i>n-i  are 
used  as  the  interpolation  functions  in  this  procedure,  they  can  be  expected  to  yield  an 
accurate  interpolation  scheme  for  band-limited  functions,  provided  that  the  matrix  to  be 
inverted  is  well-conditioned.  The  following  theorem  shows  that  if  the  interpolation  nodes 
are  chosen  to  be  quadrature  nodes  accurate  up  to  twice  the  bandwidth  of  interpolation, 
with  the  quadrature  formula  being  accurate  to  more  than  twice  as  many  digits  as  the 
interpolation  formula  is  to  be  accurate  to,  then  the  matrix  inverted  in  the  procedure  is 
close  to  being  a  scaled  version  of  an  orthogonal  matrix. 


Theorem  7.1  Suppose  the  numbers  Wi,W2,  ■  ■  ■  iWn  E  R  and  Xi,X2, ...  ,Xn  eR  are  such 
that 


<  e, 


(127) 
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for  all  a  £  [~1)  1];  for  some  c  >  0.  Let  the  matrix  A  he  given  by  the  formula 

fpo{Xi)  -tpiixi)  ...  i)n-l{xi) 

^0(2:2)  ^l{X2)  '^n-liX2) 

i’oiXn)  tpliXn)  ■■■  ^n-l{Xn) 

let  the  matrix  W  he  the  diagonal  matrix  whose  diagonal  entries  are  Wi,  W2, .  - . ,  Wn,  and 
let  the  matrix  E  =  [ejk]  be  given  by  the  formula 

E  =  I-  A*WA.  (129) 

Then 

\ejk\  <  ■  .  (130) 

Proof.  Clearly 

n 

Cjk  =  djk -'^wi'tpj.i{xi)i)k-i{xi),  -  (131) 

1=1 

where  Sij  is  the  Kronecker  delta  function.  Using  (18),  this  becomes 

e^k  =  (t—  [  e"*“'>j_i(t)  dt) 

1=1  / 

•  \T—  [  tpk-iir)  dr] 

\Ak-i  J 

=  bjk-= — - [  [  rpj-i{t)'ifk-i{r)  Ylwie~^‘^‘^  dtdr.  (132) 

•'-1  •'-i 

Using  (127),  this  becomes 

^jk  —  ^jk~=F  r  f  f  'lpk-l{x)  (133) 

Aj^lAk-i  J-l  /  rl  .  ,  .  \ 

\J  ^  ds  —  /£(t+r)j  dtdr, 

where  /e :  2, 2]  — >■  C  is  a  function  which  satisfies  the  relation 

|/e(x)l  <  e,  (134) 

for  all  X  €  [—2, 2].  Thus 

ejk  =  ^jk-^ — \ - [  [  'ipj-i{t)ipk-i{r)  [  dsdtdr 

Aj—iAk—l  •'—1  .'—1  J—1 

+=—^ -  /  [  tpj-i{t)^k-iir)  fe{t  +  T)  dtdr  (135) 

•'-1  -^-1 
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Using  (18),  this  becomes 


ejk  =  Sjk  -  ipj-iis)  ipk-iis)  ds 

+  A  "lAfe  1 1-i  Li  fS  +  'r)  dtdr. 

Due  to  the  orthonormality  of  the  functions  {•0^},  this  becomes 

^jk  =  y — \ -  [  ipk-iir)  [  fs{t  +  t)  dtdr. 

Aj-iAk-l  J-\ 

Using  the  Cauchy-Schwartz  inequality,  this  becomes 


(136) 

(137) 


(138) 

□ 


From  inspection  of  Theorem  2.5,  it  can  easily  be  seen  that  the  number  N  of  eigenval¬ 
ues  needed  for  a  bandwidth  of  2c  and  an  accuracy  of  ^  is  roughly  twice  the  number  of 
eigenvalues  needed  for  a  bandwidth  of  c  and  an  accuracy  of  e.  Thus  a  generalized  Gaus¬ 
sian  quadrature  for  a  bandwidth  2c  and  an  accuracy  has  roughly  the  same  number 
of  nodes  as  are  needed  for  interpolation  of  accuracy  e.  In  our  numerical  experiments, 
this  correspondence  was  found  to  be  much  closer  than  the  rough  bounds  in  Theorem  2.5 
indicate;  in  the  results  tabulated  in  Section  8,  the  number  of  nodes  for  an  interpolation 
formula  of  a  desired  accuracy  e  was  always  chosen  to  be  the  number  of  quadrature  nodes 
for  a  desired  accuracy  for  twice  the  band  limit  (that  number,  in  turn,  being  chosen 
as  indicated  in  Section  5);  the  correspondence  between  the  desired  accuracy  and  the 
experimentally  measured  maximum  error  can  be  seen  in  Tables  3  and  4. 

The  coefficients  Ci,C2,...,c„  produced  by  this  interpolation  procedure  (see  (126)) 
can,  of  course,  just  as  easily  be  used  for  evaluating  derivatives  or  indefinite  integrals  of 
the  interpolated  function,  as  they  can  for  computing  the  function  itself. 
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8  Numerical  Results 


The  algorithms  of  Sections  5-7  have  been  implemented  in  double  precision  (64-bit  floating 
point)  arithmetic,  with  results  shown  in  Tables  1-4.  Tables  1  and  2  show  the  perfor¬ 
mance  of  quadrature  nodes  produced  by  the  schemes  of  Sections  5  and  6,  when  used  as 
quadrature  nodes;  Tables  3  and  4  show  their  performance  when  used  as  interpolation 
nodes.  These  are  not  actually  the  same  sets  of  nodes;  even  with  the  bandwidth  c  for  in¬ 
terpolation  being  half  of  the  bandwidth  for  quadrature  (as  it  is  in  the  tables),  more  nodes 
are  needed  to  achieve  a  given  accuracy  of  interpolation  than  are  needed  to  achieve  a  given 
accuracy  of  quadrature,  as  can  be  seen  by  comparing  the  number  of  nodes  (printed  in 
the  column  labeled  n  in  each  table).  The  error  figures  in  the  tables  are  approximations 
of  the  maximum  error  of  interpolation  or  of  the  quadrature,  when  applied  to  functions 
of  the  form  cos(ax)  and  sin(orc),  with  0  <  a  <  c;  they  were  computed  by  measuring  the 
error  at  a  large  number  of  points  in  a  (for  interpolation,  in  both  a  and  x),  including  the 
extremes.  The  column  labeled  “Roots”  contains  the  errors  for  the  nodes  produced  by 
the  scheme  of  Section  6;  the  column  labeled  “Refined”  contains  the  errors  after  those 
nodes,  used  as  a  starting  point,  have  been  run  through  the  scheme  of  Section  5.  The 
variable  e  which  appears  in  the  tables  is  the  requested  accuracy,  used  to  determine  the 
number  of  nodes  in  the  ways  described  in  Sections  5  and  7. 

Also  tabulated  are  the  numbers  of  Legendre  nodes  required  to  achieve  the  same 
accuracy  e  using  polynomial  interpolation  or  quadrature  schemes.  Since  Chebyshev 
nodes  are  generally  known  to  be  superior  for  interpolation,  for  that  case  the  numbers  of 
Chebyshev  nodes  required  to  achieve  the  same  accuracy  are  also  tabulated. 

Figure  2  contains  the  maximum  norm  of  the  derivative  of  each  prolate  function  i>j{x), 
for  c  =  200  and  x  €  [-1,1],  as  a  function  of  j]  also  graphed,  for  comparison,  is  the 
maximum  norm  of  the  derivative  of  each  normalized  Legendre  polynomial  Pj{x)  over 
the  same  range;  and  graphed  below,  on  the  same  horizontal  scale,  are  the  norms  of  the 
eigenvalues  Xj.  The  graph  shows  that,  for  this  value  of  c,  computing  the  derivatives  of 
a  function  given  by  a  prolate  series  is  a  better-conditioned  operation  than  computing 
the  derivatives  of  a  function  given  by  a  Legendre  series  of  the  same  number  of  terms. 
(Obviously,  if  the  number  of  terms  can  also  be  reduced,  as  in  the  situations  of  Tables  1- 
4,  there  is  a  further  improvement  in  the  condition  number.)  The  same  general  pattern 
of  behavior  is  exhibited  for  other  values  of  c;  as  c  approaches  zero  (and  the  prolate 
functions  approach  the  Legendre  polynomials),  the  value  of  j  at  which  the  maximum 
norm  of  the  derivative  rises  sharply  also  approaches  zero  (as  is  to  be  expected,  since  for 
c  =  0  the  prolate  functions  reduce  to  Legendre  polynomials).  Finally,  Tables  5  and  6 
contain  samples  of  quadrature  weights  and  nodes. 

Remark  8.1  In  this  paper,  detailed  discussion  of  issues  encountered  in  the  implemen¬ 
tation  of  numerical  algorithms  has  been  deliberately  avoided,  as  well  as  any  discussion  of 
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CPU  time  requirements,  memory  requirements,  etc.  Thus,  we  limit  ourselves  to  observ¬ 
ing  that  all  algorithms  have  been  implemented  in  FORTRAN,  that  with  the  exception  of 
the  procedure  for  the  evaluation  of  Prolate  Spheroidal  Wave  functions  described  in  Sec¬ 
tion  4,  we  have  not  designed  or  implemented  any  new  or  original  numerical  algorithms, 
and  that  the  procedure  of  Section  4  consists  of  applying  standard  tools  of  numerical 
analysis  (diagonalization  of  a  tridiagonal  matrix)  to  the  well-known  recursion  (61).  The 
resulting  algorithm  for  the  evaluation  of  prolate  spheroidal  wave  functions  has  the  CPU 
time  requirements  proportional  to  c^,  with  a  fairly  large  proportionality  constant.  The 
procedure  of  [2],  when  applied  to  the  system  of  functions  •  •  • ,  V’zn+i  requires  order 

operations,  also  with  a  fairly  large  proportionality  constant.  On  the  other  hand,  the 
cost  of  finding  all  roots  n  of  the  function  lying  on  the  interval  [—1, 1]  is  proportional 
to  n,  and  the  proportionality  constant  is  not  large.  The  largest  c  we  have  dealt  with  in 
our  experiments  was  about  6000,  with  resulting  quadratures  having  about  1900  nodes. 
In  this  regime,  the  construction  of  the  quadrature  (both  nodes  and  weights)  took  several 
minutes  on  the  SOO-rnpgaflop  SUN  workstation;  while  there  are  fairly  obvious  ways  to 
reduce  the  cost  of  the  calculation  (both  in  terms  of  asymptotic  CPU  time  requirements 
and  in  terms  of  associated  proportionality  constants)  we  have  made  no  effort  to  do  so. 

The  following  observations  can  be  made  from  the  examples  presented  in  this  section, 
and  from  the  more  extensive  tests  performed  by  the  authors. 

1.  When  the  nodes  obtained  via  the  algorithm  of  [2]  are  used  for  the  integration  of  band- 
limited  functions,  the  resulting  quadrature  rules  are  significantly  more  accurate  than  the 
quadratures  obtained  from  the  nodes  of  appropriately  chosen  prolate  functions;  however, 
the  difference  between  the  numbers  of  nodes  required  by  the  two  approaches  to  obtain 
a  prescribed  precision  is  not  large.  When  the  nodes  obtained  via  the  two  approaches  are 
used  for  the  interpolation  (as  opposed  to  the  integration)  of  band-limited  functions,  the 
performances  of  the  two  are  virtually  identical. 

2.  For  large  c,  the  number  of  nodes  required  by  a  quadrature  rule  for  the  integration 
of  band-limited  functions  with  the  band-limit  c  is  close  to  the  dependence  on  the 
required  precision  of  integration  is  weak  (as  one  would  expect  from  Theorem  2.5  and 
subsequent  developments). 

3.  The  numbers  of  nodes  required  by  our  quadratures  rules  to  integrate  band-limited 
functions  is  roughly  7r/2  times  less  than  the  numbers  of  Gaussian  nodes;  the  numbers 
of  nodes  required  by  our  interpolation  formulae  in  order  to  interpolate  band-limited 
functions  is  roughly  7r/2  times  less  than  the  number  of  Chebychev  (or  Gaussian)  nodes. 
Again,  the  dependence  of  the  required  number  of  nodes  on  the  accuracy  requirements  is 
weak. 

4.  The  norm  of  the  differentiation  operator  based  on  our  nodes  is  of  the  order  as 


compared  to  the  norm  of  the  spectral  differentiation  operators  obtained  from  classical 
polynomial  expansions;  this  might  be  useful  in  the  design  of  spectral  (or  pseudospectral) 
techniques. 


9  Miscellaneous  Properties 

Prolate  spheroidal  wave  functions  possess  a  rich  set  of  properties,  vaguely  resembling  the 
properties  of  Bessel  functions.  This  section  establishes  some  of  those  properties.  Some 
of  the  identities  below  can  be  found  in  [20], [17], [5];  others  are  easily  derivable  from  the 
former. 


The  identity 

00 

e'cxt  =  Xj  (a;)  (t) ,  (139) 

j=0 


(see  Section  5)  has  a  number  of  consequences  which,  while  fairly  obvious,  seem  worth 
recording,  since  similar  properties  of  other  special  functions  have  often  been  found  useful. 
Differentiating  (139)  m  times  with  respect  to  x  and  n  times  with  respect  to  t  yields  the 
formula 


/  1  \  (m+n)  00 

^m^n^^cxt  ^  (140) 

for  all  x,t  £  [—1,1]-  Multiplying  (139)  by  and  integrating  with  respect  to  t, 

converts  it  into 


sin{c  ■  {x  —  u)) 
X  —  u 


00 


j=0 


(141) 


Taking  the  squared  norm  of  (139),  and  integrating  with  respect  to  x  and  t,  yields  the 
formula 


El^al'  =  4;  (142) 

j=o 

combining  this  with  (21)  yields 


^  2c 

j=o  ^ 


(143) 
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(144) 


Setting  X  =  t  =  1  converts  (139)  into 
00 

j=o 

The  identity 

Xjtjjj  (x)  =  i)j  (t)  dt  ( 145) 

(see  Section  2.5)  also  has  a  number  of  simple  but  potentially  useful  consequences.  Dif¬ 
ferentiating  it  k  times  with  respect  to  re,  we  get 

(x)  =  {ic)'^  y (t)  dt.  (146) 

We  next  consider  the  integral 

f{x)  =  f{a,x)  =  [  j - 'tpjit)  dt.  (147) 

Differentiating  (147)  with  respect  to  x,  we  have 

-^f{a,x)  =  ic  j  — ij;j{t)  dt.  (148) 

dx  J-i  t  —  a 

Multiplying  (147)  by  ica,  and  subtracting  it  from  (148),  we  obtain 

-^f{a,x)  —  icaf{a,x)  =  ic  [  ipjit)  dt  (149) 

dx  J-i 

=  icXjipj{x). 

In  other  words,  /  satisfies  the  differential  equation 

f'{x)  —  icaf{x)  =  icXjijjj{x).  (150) 

The  standard  “variation  of  parameter”  calculation  provides  the  solution  to  (150): 

f{x)  =  icXj  r  e-*“(^-‘Vj(t)  dt  /(O)  (151) 

J  0 


Introducing  the  notation 


V  = 


1  d 
ic  dx 


(152) 
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(i.e.  V  is  the  product  of  multiplication  by  1/ic  and  differentiation),  we  rewrite  (146)  as 


=  ^  r  (153) 

Aj  J-1 

for  an  arbitrary  polynomial  P  (with  real  or  complex  coefficients), 

=  i  £  P(i)  ■(«)*.  (154) 

By  the  same  token,  the  function  ^  defined  by  the  formula 

(155) 

satisfies  the  differential  equation 

(156) 


The  following  lemma  provides  a  recursion  connecting  the  values  of  the  A;-th  derivative 
of  the  function  with  its  derivatives  of  orders  k  —  1,  A;  —  2,  A:  —  3,  A;  —  4. 

Lemma  9.1  For  any  positive  real  c,  integer  m  >  0,  and  x  €  (— oo,  +oo), 


(1  -  2;^)  -  2  (A:  +  1)  a;  ^^+^^(a:) 

+  (Xm  -  A:  (A:  +  1)  -  x^) 

-  2  (?  k  X  ip^~^\x)  —  (?k{k  -1)  '(pl^~‘^\x)  =  0  (157) 

for  all  k  >2.  Furthermore, 

(1  -  ^m(^)  -  (a;)  4-  (Xm  -  2  -  x^)  'if'„,{x) 

-  2  (?  X 'lpm{x)  =  Q  .  (158) 

In  particular, 

-2{k  +  l)  +  (x™  -  M*:  +  1)  -  c")  V-g*  (1) 

-  2  cH  </>£-')  (1)  -<fk{k-  1)V'S-*>(1)  =  0  (159) 

for  all  k>2,  and 

-  2  V'm(l)  +  (Xm  -  C^)  V’m(l)  =  0  ,  (160) 

-  4 C(l)  +  (Xm  -  2  -  c^)  -  2  c2  ^^(1)  =  0 .  (161) 
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Furthermore,  for  all  integer  m  >  0  and  k>2, 

V'r=>(0)  +  (x™-MA:  +  l))V'S>(0) 


(162) 

For  all  odd  m, 

C(0)  +  (Xm-2)C(0)=0. 

(163) 

and  for  all  even  m, 

^m(O)  +  Xm'i/''m(0)  =0. 

(164) 

Finally,  for  all  integer  m  >  0,  /c  >  0, 

(165) 

'/>Sii(o)=o, 

(166) 

V-S*')(o)  =  o. 

(167) 

Proof.  All  of  the  identities  (157)  -  (164),  (166),  (167),  are  immediately  obtained  by 
repeated  differentiation  of  (24). 

In  order  to  prove  (165),  we  assume  that 


M)  =  0  (168) 

for  some  integer  m  >  0.  and  observe  that  the  combination  of  (168)  with  (159),  (160),  (161) 
implies  that 

V’S>(1)  =  0  (169) 

for  all  A:  =  0, 1, 2, ... .  Due  to  the  analyticity  of  t^m(^)  in  the  complex  plane,  this  would 
imply  that 


'(pm{x)  =  0 

for  all  a;  € 

The  following  is  an  immediate  consequence  of  the  identity  (160)  of  Lemma  9.1. 
Corollary  9.2  For  all  integer  m,  n  >  0, 


(170) 

□ 


^m(l)  •  '0n(l)  -  •0n(l)  '  ^/>m(l)  =  (Xn  “  Xm)  •  '0n(l)  '  V’m(l)  ,  (171) 


where  Xm,Xn  ^  R  are  as  defined  in  Theorem  2.6. 
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Theorem  3.1,  in  Section  3.1,  gives  formulae  for  the  entries  of  matrices  for  differen¬ 
tiation  of  prolate  series  and  for  multiplication  of  prolate  series  by  x.  Matrices  for  any 
combination  of  differentiation  and  of  multiplication  by  a  polynomial  can  obviously  be 
constructed  from  these  two  matrices;  for  instance,  calling  the  differentiation  matrix  D, 
and  the  multiplication-by-x  matrix  X,  the  matrix  for  taking  the  second  derivative  of  a 
prolate  series,  then  multiplying  it  by  5  —  is  equal  to  (5/  —  X'^)D^. 

In  many  cases,  however,  there  are  simpler  formulae  for  the  entries  of  such  matrices, 
that  is,  for  inner  products  of  ipj{x)  with  its  derivatives  and  with  polynomials.  The  follow¬ 
ing  theorems  establish  several  such  formulae,  as  well  as  a  few  formulae  for  inner  products 
which  do  not  involve  'ipjix)  itself  but  only  its  derivatives.  We  start  with  Theorem  3.1, 
restated  here  for  consistency. 

Theorem  9.3  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  If  m  =  n  (mod  2),  then 

j  ^  'tf'^{x)  Xpmix)  =  J  'ipmix)  dx  =  0.  (172) 

Ifm^n  (mod 2),  then 

[  Ip'nix) 'Ipmix)  dx  =  .2^X2  ^m(l)  V^n(l)  ,  (173) 

•/-I  +  K 

[  X'4)n{^)'(pm{x)  dx  =  ^  Tprnil)  ^n{l)  ■  (174) 

J-1  ^c  -h 

Theorem  9.4  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  Ifm^n  (mod 2),  then 

X  i)'^{x)  -ipTnix)  dx  =  0.  (175) 

Ifm  —  n  (mod  2),  then 

[  x-ip'^ix)  dx  =  (2'0m(l)  ^n(l)  “  ^mn)  •  (176) 

Proof.  Identity  (175)  is  obvious  since  the  functions  ijjj  are  alternately  even  and  odd  (see 
Theorem  2.4).  In  order  to  prove  (176),  we  consider  the  integral 

j  ^X^l)!,,{x)i;rn{x)  dx 

"T  /  ^([  dt)  i}m{x)  dx 

An  J-1  \J-1  J  X 
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(177) 


=  ^  J  ^xil>m{x)  (^J  dt^  dx 

—  ^  X  tprni^)  dx^  1pn{t)  dt 

=  ^  j_^t'>P'mif)'^n{t)  dt. 

In  other  words, 

fl  rl 

/  X  'tpm{x)  dx=~  X  'Ip'^ix)  ^n(x)  dx  . 

J-1  Aji  J-l 

On  the  other  hand,  integrating  the  left  side  of  (177)  by  parts,  we  obtain 
x'ip'.^{x)ipm{x)  dx 

=  2i)rn{l)^n{l)  -  {'ljjn{x)  ^p!^{x)  X  +  tjjn{x)  ^mix))  dx 

=  2lprn{l)'lpnil)  -  X 'ipn{x)  i^'mix)  dx  -  5ran  ■ 

Combining  (177)  and  (178),  we  have 

xip'^{x)ipn{x)  dx 

/\fi  y  —1 

=  2^^(l)  1pn{l)  ~  ^  ^m(^)  i’nix)  dx  -  5mn  , 

from  which  (176)  follows  directly.  □ 


Theorem  9.5  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  Ifm^n  (mod 2),  then 


j  ^  i’mix)  dx  =  0  . 


(178) 


Ifm  =  n  (mod  2)  and  m  ^  n,  then 

pI 

j_  fpnix)  dx 

2\ 

=  3-f^«(l)^7.(l)-C(l)V'n(l)) 


(179) 
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(180) 


'  "  ,  {Xn  -  Xm)  ^(1) 


4An 


'0n(l) 


An  Am 

where  XmjXn  ^  R,  are  as  defined  in  Theorem  2.6. 


Proof.  Clearly  (178)  is  true,  since  the  functions  ■0^  are  alternately  even  and  odd.  In 
order  to  prove  (179)  and  (180),  supposing  that  m  =  n  (mod  2)  and  m  ^  n,  we  consider 
the  integral 

rl 

~  Y  '  il-l  '^rn{x)  dx 

=  J  ^  fi}m{x)x^-  ^  dt^  dx 

=  (^J  ^  'ipm{x)  dx'^  1pn{t)t^  dt 

=  ^  f  t2^„(t)^"  (t)  dt, 

/\j^  J  —1 

which  is  summarized  as 

[  X^  1p!^{x)  iJrn{x)  dx  =  ^  [  X^  Tp'^x)  i}n{x)  dx  .  (181) 

J —  1l  Ayj  J —  1 

On  the  other  hand,  integrating  the  left  side  of  (181)  by  parts,  we  have 

rl 

X^  i}!^{x)  tprn{x)  dx 

=  2V'Ul)^m(l)  -  (i}'Y^)x^  +  2xi)rn{x))  dx 

=  2^p'^{l)'lpm{l) -2  i>'^{x)^prnix)x  dx 

-  j_^i^'n{^)i>'m{x)x^  dx.  (182) 

Due  to  Theorem  9.4  and  the  fact  that  n,  we  immediately  rewrite  (182)  as 

rl 

X^ 'lp'^{x) 'tprnix)  dx 

2\ 

=  2^;(l)^m(l)  -  '  -fx  2V>n(l)^m(l) 

-  dx,  (183) 
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which  we  rewrite  as 

x'^ 'ip'^ix) 'tp'^ix)  dx 

4  A 

=  2  i;rn{l)  -  i)nil)  ^m(l) 

Ajti  “t  A71 

-  x'^  1p!^{x)  i^rn{x)  dx  .  (184) 

Swapping  m  with  n,  we  convert  (184)  into 

x^'ip'^{x)ip'^{x)  dx 

=  2  C(l)  ^n(l)  -  --j"  ■  -  'IpniX)  ^m(l) 

-  x'^  ip!^{x)  i^nix)  dx.  (185) 

Combining  (184)  and  (185),  we  obtain 

/•l  «  4 

/  <(x)  tprnix)  dx  -  2  ^m(l)  +  .  ,  T"  ^n(l)  V'm(l) 

7  —  1  Atti  "T”  Ayj 

rl  AX 

=  /  x2t/>"(x)t/)„(x)  dx-2V’(„(l)V’n(l)+ 

7  —  1  Am  ‘  Ati 

which  is  obviously  equivalent  to 

X^  tp'^ix) 'tprnix)  dx 

=  j_^  X^  ij'^ix)  '4)n{x)  dx  +  2  V’m(l)  “  ^n(l)) 

+  4^^^n(l)V’m(l). 

An  “r  Ayn 

Finally,  combining  (181)  with  (187),  we  have 
^  (x)^„(x)  dx 

An  V— 1 

=  1?  ’l’a(x)  dx  +  2  «(1)  tf„(l)  -  ViiCl)  </'«(!)) 

+^r^Wi)«i).  (187) 

An  H“  Ayn 
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which  is  easily  rewritten  as 

f 


1 


dx 


=  2  (*(l)V>„(l)-C(l)^'’.(l)) 

+  4  ^”*  IpnCl)  t/'m(l)  , 


or 


2A„ 


dx 


Am  Ayj 
4  An 


An  4“  Aj^ 


(V'n(l)^m(l)-'0m(l)^n(l)) 


(188) 


We  finally  rewrite  (188)  as  (180)  using  Corollary  9.2.  □ 

The  following  theorem  is  an  immediate  consequence  of  combining  the  preceding  theorem 
with  equation  (184)  from  its  proof. 

Theorem  9.6  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  If  m  ^  n  (mod  2),  then 


j  j  v;(x)  i;'^{x)  dx  =  0 . 
If  m  =  n  (mod  2)  and  m  ^  n, 

r\ 

j  x‘  v'„(x)ii'^(x)  dx 


(189) 


=  2t;4(l)V-„(l)  + 


2A„ 


=  2<(1),).„(1)  + 


Am  An 

m 


(C(l)^n(l)  -'0n(l)V'm(l)) 


2A, 


An  An 


(V'n(l)^m(l)  -  V’m(l)V’n(l)) 


=  ti'’m(l)^n(l) 


'n  ^'m 

AmXm  AnXn 
Anx  An 


-  0“ 


(190) 

(191) 

(192) 


Theorem  9.7  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m  and  n  are 
non-negative.  Ifm^n  (mod  2),  then 


-ifnix)  t^m (x)  dx  =  X^  Ipnix)  i’mix)  dx  =  0 


(193) 
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If  m  =  n  (mod  2)  and  m^n,  then 
[  i^nix) 'ip'^ix)  dx 

J  —  X 

9 


J  ^X^1pn{x)'lpm{x)  dx 

(V'n(l)^m{l)  -^n(l)C(l)) 

iXn  Xm)  '0m(l)  '0n(l)) 


2 

1  AttjAtj 

'?  A2  -  A2 


w/iere  Xm,Xn  €  R  are  as  defined  in  Theorem  2.6. 

Proof.  Identity  (193)  is  obvious,  since  the  functions  'ifj  are  alternately  even 
In  order  to  prove  (194)-(197),  we  start  with  the  expression 


e - ■0„(t)  dt. 

Taking  the  inner  product  of  (198)  with  iprnix),  we  have 
A. 


which  we  summarize  as 


'fp'n{x)'ipm{x)  dx 

=  j  ^  (^J  ^  dt^  'tprnix)  dx 

=  (^j  ^  '4)rn{x)e'‘'^*  dx'^  dt 

rl 

>^mj  ^  i^^n(i)  V’m(i)  dt , 


/I  1  A 

'tpnix)  tpmix)  dx  = - T  /  t/;"(x)  lfm{x)  dx  . 

-1  Xm  j-l 


•1  Ajjj 

Swapping  n,  m,  we  rewrite  (199)  in  the  form  of 


J^^X^ 'ifnix)  1pmix)dx 

=  -4^  /,  '>Pm{x)'lfn{x)  dx. 

C  Aji  */  — 1 


(194) 

(195) 

(196) 

(197) 

and  odd. 

(198) 


(199) 

(200) 
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(201) 


Combining  (199)  and  (200),  we  get 

tpmioo)  dx  =  ^  ip'^{x)  'tpni^)  dx . 

On  the  other  hand,  integrating  the  left  side  of  (201)  by  parts,  we  have 

J_^'>Pn{^)  i’mix)  dx 

=  V'm(l)  -  '0n(-l)  ^m(-l)  ”  dx 

=  2  t/'n(l)  1pm{l)  -  (V'n(l)  “  V’n(-l)  ^(“1)) 

+  'tpn{x)ip'^{x)  dx.  (202) 

We  rewrite  (202)  in  the  form  of 

ipl{x)  ipm{x)  dx 

=  2  Xprnil)  -  '0n(l)  ^n{x)  1p'^{x)  dx  . 

We  combine  (201)  and  (203)  and  get 

-  ')  /j  'Pn(x)i>!n(x)  dx 

=  2  (*(l)«l)-'/'»(l)  *(!)).  (203) 

Since  m  ^  n,  we  easily  rewrite  (203)  as  (194).  We  obtain  expression  (196)  by  combin¬ 
ing  (200)  and  (194).  The  identities  (195),  (197)  follow  from  (194),  (196)  immediately 
due  to  Corollary  9.2.  □ 


Theorem  9.8  Suppose  that  c  is  real  and  positive,  and  that  the  integers  m 
non-negative.  Let 

and  n  are 

ry 

'J'n(y)  =  /  ‘tpn{x)dx. 

Jo 

(204) 

If  n  is  odd  and  m  is  even,  then 

j  ^  j’Pn(t)^m(t)  dt 

(205) 

=  ic^^^»„(l)»„(l) 

'^n 

(206) 

XI /-I 

(207) 
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(208) 


If  m  =  n  (mod  2),  then 

dt  =  0. 

Proof.  We  start  with  the  identity 

An  ^n  (a:)  =  J  ^  V’n  (i)  dt .  (209) 


Integrating  (209)  with  respect  to  x,  we  have 
ry 

An  /  'tpn{x)dx 
Jo 


=  ^  (^J  ipn{t)  dt^  dx 

(210) 

=  J  fpn{t)  j  da;dt 

(211) 

1  /•!  1  1 

=  -  -Mt)e^^'dt-- 

ZC  7-1  t  ZC  J 

[  ^i’n{t)dt, 

—  1  t 

(212) 

which  we  summarize  as 

An  ^n(y)  =  —  /  7  fpriit)  dt 

ZC7-1  t 

1  /■!  1 

-—  /  -'tpn{t)dt. 
zc  7-1  t 

(213) 

Taking  the  inner  product  of  (213)  and  ipmiy),  we  obtain 

=  h  /-,  ■  (£  1  *) 

~icLi  ’  (/i  7 

-— [  ^^n{t)dt-f  ipm{y)  dy  (215) 

2C»/— Iv  y  —  1 

=  ^/  dt 

%  C  j  •“  1  V 

/  ^'tpn{t)dt-f  ipm{y)dy,  (216) 
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which  we  summarize  as 


J  Ipnit)  Ipmit)  dt 


=  iC^  [  '^n{t)i^m{t)  dt 

J —  X 

1  /•!  1  /•! 

+  7—  /  7^n(i)  dt  -  /  dy 

^77T  7— 1  t  7  —  1 

(217) 

Exchanging  m  with  n,  we  convert  (217)  into 

j  ^  ^  -IpTnit)  Xpn{t)  dt 

=  dt 

Afi  7—1 

1  /•!  1  /*! 

+  ^  y  1  -  xpm{t)  dt-  j  ^  1pn{y)  dy  , 

(218) 

and  combining  (217),  (218),  we  get 

f  dt-^ic  [  ^m{t)'ipn{t)  dt 

■^m  J —  X  Aji  J —  1 

=  ^  [  dt  -  [  ^nit)  dt 

Afl  J —  1  Z  J —  I 

[  7^n(i)  dt  -  [  iprn{t)  dt .  (219) 

Ajji  •/  —  I  Z  J —’1 

Suppose  that  m  is  even  and  n  is  odd;  then  the  first  product  in  the  right  hand  side  of  (219) 
is  zero,  so 

^ic  [  ^n{t)  iXmit)  dt-^ic  [  ^rnit)  Xpn{t)  dt 

Am  J  Ayi  7—1 

=  f  ];xpn{t)dt-[  'tpmit)dt,  (220) 

Am  7  — 1  b  7  —  1 

which  is  equivalent  to 

/  ^n(i)  xjjmit)  dt 

J  —1 

=  ^  /  '^mit)i)n{t)  dt 

^i)n{t)dt-f  1pra{t)dt,  (221) 

An  ^C  J-1  t  J-1 
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or 


<ifrn{t)tpn{t)  dt 
=  ^  /  dt 

+  lfPn{t)dt-[  1prn{t)dt.  (222) 

/\j^  "t  C  «/”l  C  «' “1 

On  the  other  hand,  integrating  the  left  side  of  (222)  by  parts,  we  obtain 
'^m{t)ipn{t)  dt 

=  ^'n(l)  ^^,„(1)  -  ^„(-l)  «'m(-l)  -  ^n(t)  ^m(i)  dt .  (223) 


Since  the  product  ^^n(^)  is  an  odd  function  when  m^n  (mod  2),  we  rewrite  (223) 


as 


'lpn{t)  dt 

=  dt. 


The  combination  of  (222)  and  (224)  implies  that 

A?  /■! 


[  '^n{t)tpm{t)  dt+  ^  [  '^n{t)lpm(t)  dt 

J-1  J-1 

=  2  ^'„(1)  ^rn{l)  7  ’ 


or 


A^  +  A^  /•! 


A2 

m 


which  is  equivalent  to 

'■i 


2X1 


A^  +  A^ 


«^„(1)^,„(1) 


^  /  7^n(t)  dt  -  [  ^rnit)  dt . 
t  0  J  —~X  L  «/'”l 


A^  + A^  ic 


(224) 


(225) 


’l'n(0^m(t)  dt 

[  ^1pn{t)dt-[  7prn{t)  dt ,  (226) 

2  C  J-1  t  J-1 


(227) 
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Finally,  combining  (217)  and  (227),  we  have 


2  Am  A. 


=  IC 


Am 


dt 

^n(l)^m(l) 


A2+A2 


J  i’nit)  dt  ■  'Iprnit)  dt . 


(228) 


Equation  (208)  is  easily  proven  since  the  product  |  'ijjm{x)  ipn{x)  is  an  odd  function  when¬ 
ever  m  =  n  (mod  2).  □ 


The  above  theorems  do  not  use  much  of  the  detailed  structure  of  the  integral  operators 
of  which  the  functions  {ipj}  are  eigenfunctions.  Thus  many  of  them  generalize  easily  to 
the  case  of  an  operator  L  :  L^[0, 1]  — >■  L'^[0, 1]  defined  via  the  formula 

L{v)(i)  =  [  K{xt)'tp{t)  dt,  (229) 

Jo 

for  some  function  A'  :  [0. 1]  C;  the  following  theorem  is  an  example  of  this. 

Theorem  9.9  Let  A1.A2  be  two  eigenvalues  of  the  operator  L  defined  by  (229),  that  is, 


[  Kiitjvift)  dt  =  Xi'ipi{x), 

Jo 

(230) 

[  K[xt)v2{t)  dt  =  A2^2(a:). 

Jo 

(231) 

Then 

A2  xil)[{x)il)2{x)  dx 

f  X  t/>2 (x) '01  (x)  dx 

Jo 

(232) 

provided  that  neither  Aj  nor  the  denominator  of  the  right  hand  side  of  (232)  is  zero. 
Proof.  Differentiating  (230),  (231)  with  respect  to  x,  we  get 


[  t  K' (xt) 'ipi{t)  dt  =  Ai0((x), 

Jo 

(233) 

[  t  K' (xt)  ip2{t)  dt  =  A202(a;). 

Jo 

(234) 

Multiplying  (233)  by  x02(x),  we  have 

Ai  x0'i(x)  02(x)  =  x02(x)  j  tK'{xt)ipi{t)  dt. 

(235) 
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Integrating  on  the  interval  [0, 1],  we  obtain 

Ai  [  x'tp[{x)ip2{x:)  dx  =  f  x'iIj2{x)  f  tK'{xt)'ipi{t)  dtdx  (236) 
Vo  Jo  Jo 

=  t'tpi{t)  j  X  K'{xt)tp2{x)  dxdt.  (237) 

Renaming  the  variables  of  integration  on  the  right  hand  side  from  x  to  t  and  vice  versa, 
we  get 

Ai  [  X 'ip[{x)  ip2{x)  dx  =  [  xi)i{x)  [  tK'{xt)'ip2it)  dtdx.  (238) 

Jo  Jo  Jo 

Substituting  (234)  into  (238),  we  obtain 

Ai  f  X 'ip[{x) '4)2{x)  dx  =  X2  f  xxpi{x)ip2{^)  dx,  (239) 

•v  0  Jo 

from  which  (232)  follows  immediately,  as  does  its  caveat.  □ 

The  following  theorem  establishes  the  relation  between  the  norm  of  each  function  tpj 
on  [-1, 1]  (which  in  this  paper  is  taken  to  be  one),  and  its  norm  on  (—00, 00). 

Theorem  9.10  Suppose  that  c  is  real  and  positive,  and  that  the  integer  n  is  non¬ 
negative.  Then 


/OO  1 

V’n(^)  dx=—. 

-00  fJ/fi 


(240) 


where  is  given  by  (21). 

Proof. 


sin(c  •  {x  —  t)) 


X  —  t 


Pn  V-1  \7r  J-00  X  —  t  j 

=  y  [  ylit)  dt 

fJn  •'—1 


dt 


Pn 

Pn 


□ 


The  following  theorem  extends  Theorem  (9.10)  to  any  band-limited  function  with 
band  limit  c. 
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Theorem  9.11  Suppose  that  c  is  real  and  positive,  that  the  integer  n  is  non-negative, 
and  that  /  :  R  — >■  C  is  a  band-limited  function  with  band  limit  c.  Then 


Proof. 


/oo  1 

tpnix)  f{x)  dx  =  —  f{x)  dx. 

•OO  jin  ‘'“I 


/oo 

i}n{x)f{x)  dx 

-00 

7-00  yKPn  7-1  x  —  t  J 

7-1  \7r  7-00 


x  —  t 


—  [  'lpn{t)f{t)  dt. 

fJ>n  7  —  1 


(241) 


□ 


Theorem  9.12  Suppose  that  c  is  real  and  positive,  and  that  the  integer  n  is  non¬ 
negative.  Then 


/OO 

-00 


—'ifrnix),  if  -  I  <X  <1, 


0, 


if  X  >  1  or  X  <  — 1. 


(242) 


Proof.  Since  ifrn  is  an  eigenfunction  of  the  operator  Qc  defined  in  (19),  and  pm  is  the 
corresponding  eigenvalue, 


Thus 


Mm'0m(^)  ~  / 

TT  J-l 

/OO 

dt 

-00 


1  sin{c  ■  {x  —  u)) 


X  —  u 


du. 


_  1  f°° 

U>m  — ' 


/-/“M  [If 


1  /■!  sin{c  ■  [x  —  u)) 


l^m  oo 

_L 

Pm 


X  —  U 


X  —  u 


(243) 


1prn[u)  du  j  dt 

(244) 

/ 

^e^^^dt]  du 

(245) 
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Since  the  innermost  integral  is  the  orthogonal  projection  operator  onto  the  space  of 
functions  of  band  limit  c  on  (—00,00),  applied  to  the  function  it  follows  that: 


/: 


dt 


_  J  l^m 


0, 


if  —  1  <  x  <  1, 
if  X  >  1  or  X  <  — 1 


1  /•! 

=  -  /  V'm(M) 

\  —  f  ipmiu)  du,  if  -  1  <  X  <  1, 

J  Mm  J~“l 

[0,  if  X  >  1  or  X  <  — 1, 


du 


(246) 

(247) 


from  which  (242)  follows  immediately. 


□ 


The  following  five  theorems  establish  formulae  for  the  derivatives  of  prolate  functions 
and  their  associated  eigenvalues  with  respect  to  c. 

Theorem  9.13  For  all  positive  real  c  and  non-negative  integer  m, 


,  2^5,(1)-1 


dc 


—  At, 


2c 


(248) 


Proof.  We  start  with 


ATn^m(2;)  =  J  ^  1pm{t)  dt . 


(249) 


Differentiating  (249)  with  respect  to  c,  we  obtain 

^Ajtj  /  \  I  \ 

Xr, 


dc 


-  L 


dc 

ixte^'^^'ipmit)  dt-\- 


d-ipruji) 

dc 


dt . 


(250) 


Multiplying  by  on  both  sides  of  (250),  and  integrating  on  the  interval  [—1, 1],  we 

get 


'ipm{x)  +  Xm - 5 - 


=  J  'tpm{x)  J  ixte'''^^ 'tprnit)  dt  dx 


(251) 
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which  we  rewrite  as 


dXm 

dc 


+  A,- 


£ 


dtpmjx) 

dc 


dx 


=  J  J  X 'lpra{x)  dx  dt 


dx  dt 


f  1 

—  Xm  I  - 

J-1  I C 


1  d'lprni't) 


dc 


dt 

Iprnit)  dt, 


dt 


which  we  summarize  as 
dXm 


dc 


On  the  other  hand,  integrating  the  right-hand  side  of  (254)  by  parts,  we  have 

which  we  rewrite  as 


(252) 


(253) 


(254) 


(255) 


t'lpmit) 


Finally,  substituting  (256)  into  (254),  we  get 
dXm  ,  2^^(1)-1 


dc 


An 


2c 


(256) 

(257) 

□ 


Theorem  9.14  For  any  positive  real  c  and  non-negative  integer  m, 


dpin 

dc 


2 

~  '0m  (1)  • 


(258) 
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(259) 


Proof.  We  start  with  the  identity 


A^m  —  '^m  '^m- 

TT 

Differentiating  (259)  with  respect  to  c,  we  get 


9fJ>n 


TT 


\r] 


dX, 


+  Xn 


dXr. 


r\  \  '  'tn  rv  '  '  'fit,  r\ 

OC  TT  \  OC  OC 

Substituting  Lemma  9.13  into  (260),  we  get 
dfx. 


"H  Am  Ar; 


TT 


dc 


2c  ,  2ii,{l)-l  ,  2-  , 

•  Z  Am  Am  ^  "r  Am  Am 

TT  2C  TT 


—  2 


2V'^(1) 


1  1 
“f"  (Xm 
C 

1 


2c 

2  1  ... 

~  A^m  '0m(^)  A^n^  H  l^m 
c  c  c 

2 

—  //m  V'm(l)  ■ 


(260) 


(261) 


(262) 

□ 


The  following  theorem  immediately  follows  from  Theorems  9.13  and  9.14. 

Theorem  9.15  For  all  positive  real  c  and  non-negative  integer  m,  n, 

=  ■  (264) 

\  J  l^n  ^ 

Theorem  9.16  Suppose  that  c  is  real  and  positive,  and  the  integers  m,n  are  non¬ 
negative.  If  n,  then 

[  ^m{t)  ^(t)  dt  =  --  ■  ^m(l)  ^n(l)  ■  (265) 

J-1  dc  c  A^  -  A^ 

If  m  =  n,  then 

I-i 
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Proof.  Since  the  norm  of  ipn  on  [—1, 1]  remains  constant  as  c  varies,  must  be  orthog¬ 
onal  on  [—1, 1]  to  its  own  derivative  with  respect  to  c,  which  immediately  yields  (266). 
To  establish  (265),  we  start  with  the  identity 


Xn'lpnix)  =  j  ^  dt  . 


(267) 


Differentiating  (267)  with  respect  to  c,  we  get 
■Ipnixj+X, 


dc 


"  dc 

=  j  +  e 


icxt  „/.  /'4\  I  ^icxt  dignity 


dc 


dt . 


(268) 


Multiplying  both  sides  of  (268)  by  and  integrating  with  respect  to  x,  we  have 


A, 


.  J  ^  '4^m{x) 


dlpnjx) 

dc 


dx 


An 

C 


X  Ip'nix)  dx  -h  Xm  'tpmit)  ,  (269) 


which,  using  (176),  we  rewrite  as 


(A„  -  A„i)  ijjmit)  dt 


An  Ar 


(2V^^(l)V^n(l)  dmri)  ■ 


C  X-m  4"  An 

Assuming  that  mj^n,  and  dividing  by  An  —  A^,  we  then  get  (265). 


(270) 

□ 


Theorem  9.17  Suppose  that  c  is  real  and  positive,  and  the  integer  m  is  non-negative. 
Then 

^  =  2cy^^2:2^^(x).  (271) 

Proof.  Due  to  Theorem  2.6, 

(1  -  x‘^)'i(j!^{x)  -  2x'lp'^{x)  +  (Xm  -  C^x^)  'lpra{x)  =  0.  (272) 

Making  the  infinitesimal  changes  c  —  c  + h,  Xm  —  Xm  +  and  i^rn{x)  =  i’mix)  +  <5(a;), 
this  becomes 


(1  -  x^)  ■  (^"  (x)  5"{x))  -  2x  ■  (ip'^ix)  -f  (5'(x)) 

+  (Xm  +  £  -  (c  +  h)  V)  •  {^pm{x)  +  5(x))  =  0.  (273) 
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Expanding  each  term,  discarding  infinitesimals  of  the  second  order  or  greater  (that  is, 
products  of  two  or  more  of  the  quantities  h,  e,  and  5(rr)),  and  subtracting  (272),  we  get 


(1  —  x^)  5"{x)  —  2x5' {x)  +  (xm  —  5{x)  +  (e  —  2chx^)ipjn{x)  =  0.  (274) 

Let  the  self-adjoint  differential  operator  L  be  defined  by  the  formula 

=  (1  -  -  ‘^xf'{x)  +  (xm  -  cV)/(a;).  (275) 

Then,  multiplying  (274)  by  tj)m{x)/h  and  integrating  on  [-1, 1],  we  get 

I-i  ^  I  “  =  0-  (276) 

Now  I  =  In  addition,  since  L  is  self-adjoint, 

L{^^){x)  'Ipmix)  dx  =  L{‘^m){x)  dx.  (277) 

But  due  to  (272),  L{'tpm){x)  =  0  for  all  x  G  [—1,1],  so  the  integral  (277)  is  zero. 
Thus  (276)  becomes 

^  =  2cy^^a;2v^^(a:).  (278) 

□ 


10  Generalizations  and  Conclusions 

In  this  paper,  we  design  quadrature  rules  for  band-limited  functions,  based  on  the  prop¬ 
erties  of  Prolate  Spheroidal  Wave  Functions  (PSWFs),  and  the  connections  of  the  latter 
with  certain  fundamental  integral  operators  (see  (17),  (19)  in  Section  2.5).  The  quadra¬ 
tures  are  a  surprisingly  close  analogue  for  band-limited  functions  of  Gaussian  quadratures 
for  polynomials,  in  that  they  have  positive  weights,  are  optimal  in  the  appropriately  de¬ 
fined  sense,  and  their  nodes,  when  used  for  approximation  (as  opposed  to  integration), 
result  in  extremely  efficient  interpolation  formulae.  Thus,  Sections  5-7  of  this  paper  can 
be  viewed  as  reproducing  for  band-limited  functions  much  of  the  standard  polynomial- 
based  approximation  theory  (for  which  see,  for  example,  [24]).  Generally,  there  is  a 
striking  analogy  between  the  band-limited  functions  and  polynomials. 

Obviously,  there  are  certain  differences  between  the  resulting  apparatus  and  the  stan¬ 
dard  numerical  analysis.  To  start  with,  where  the  classical  techniques  are  optimal  for 
polynomials,  the  approach  of  this  paper  is  optimal  for  band-limited  functions;  whenever 
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the  functions  to  be  dealt  with  are  naturally  represented  by  trigonometric  expansions  on 
finite  intervals,  our  quadrature  and  interpolation  formulae  tend  to  be  more  efficient  than 
those  based  on  the  polynomials.  When  the  functions  to  be  dealt  with  are  naturally  rep¬ 
resented  by  polynomials,  the  classical  approach  is  more  efficient;  however,  many  physical 
phenomena  involve  bend-limited  functions,  and  very  few  involve  polynomials. 

Qualitatively,  the  quadrature  (and  interpolation)  nodes  obtained  in  this  paper  behave 
like  a  compromise  between  the  Gaussian  nodes  and  the  equispaced  ones:  near  the  middle 
of  the  interval,  they  are  very  nearly  equispaced,  and  near  the  ends,  they  concentrate 
somewhat,  but  much  less  than  the  Gaussian  (or  Chebychev)  nodes  do.  For  large  c,  the 
distance  between  nodes  near  the  ends  of  the  interval  is  of  the  order  with  the  total 
number  of  nodes  close  to  ;§.  In  contrast,  the  distance  between  the  Gaussian  nodes  near 
the  ends  of  the  interval  is  of  the  order  with  n  the  total  number  of  nodes.  A  closely 
related  phenomenon  is  the  reduced  norm  of  the  differentiation  operator  based  on  the 
prolate  expansions:  for  an  n— point  differentiation  formula,  the  norm  is  of  the  order 
as  opposed  to  v?  for  polynomial-based  spectral  differentiation.  Thus,  PSWFs  are  likely 
to  be  a  better  tool  for  the  design  of  spectral  and  pseudo-spectral  techniques  than  the 
orthogonal  polynomials  and  related  functions. 

Much  of  the  analytical  apparatus  we  use  was  developed  more  than  30  years  ago 
(see  [20]- [21],  [17],  [18]);  the  fundamental  importance  of  these  results  in  certain  areas  of 
electrical  engineering  and  physics  has  also  been  understood  for  a  long  time.  However, 
there  appears  to  have  been  no  prior  attempt  made  to  view  band-limited  functions  as  a 
source  of  numerical  algorithms.  Generally,  there  is  a  fairly  limited  amount  of  information 
in  the  literature  about  the  PSWFs,  especially  when  compared  to  the  wealth  of  facts  on 
many  other  special  functions.  Section  9  of  this  paper  is  an  attempt  to  remedy  this 
situation  to  a  small  degree. 

The  apparatus  built  in  this  paper  is  a  strictly  one-dimensional  one.  Obviously,  one 
can  construct  discretizations  of  rectangles,  cubes,  etc.  by  using  direct  products  of  one¬ 
dimensional  grids;  the  resulting  numerical  algorithms  are  satisfactory  but  not  optimal. 
Furthermore,  representation  of  band-limited  functions  on  regions  in  higher  dimensions 
is  of  both  theoretical  and  engineering  interest.  Obvious  applications  include  seismic 
data  collection  and  processing,  antenna  theory,  NMR  imaging,  and  many  others.  When 
the  region  of  interest  is  a  sphere,  most  of  the  necessary  analytical  apparatus  can  be 
found  in  [21].  At  the  present  time,  we  have  constructed  and  implemented  somewhat 
rudimentary  versions  of  the  relevant  numerical  algorithms;  we  are  conducting  numerical 
experiments  with  these,  and  will  report  the  results  at  a  later  date.  A  much  more  difficult 
set  of  questions  is  presented  by  the  structure  of  band-limited  functions  on  more  general 
regions. 
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Table  1:  Quadrature  performance  for  varying  band  limits,  for  e  =  10  ^ 


c 

Maximum  Errors 

■^pol 

Roots 

Refined 

10.0 

9 

0.96E-05 

0.51E-07 

13 

20.0 

13 

0.17E-04 

0.94E-07 

19 

30.0 

17 

0.12E-04 

25 

40.0 

20 

31 

50.0 

24 

0.35E-05 

0.83E-07 

37 

60.0 

27 

0.25E-04 

0.27E-06 

43 

70.0 

31 

O.llE-04 

48 

80.0 

34 

0.48E-05 

0.17E-06 

54 

90.0 

38 

0.21E-05 

59 

100.0 

41 

0.12E-04 

0.91E-07 

65 

200.0 

74 

0.24E-05 

118 

300.0 

106 

0.32E-05 

0.21E-06 

171 

400.0 

139 

0.52E-05 

0.62E-07 

223 

500.0 

171 

0.88E-07 

275 

600.0 

203 

0.58E-05 

O.llE-06 

326 

700.0 

235 

0.12E-06 

377 

800.0 

267 

0.55E-05 

0.13E-06 

428 

900.0 

299 

0.14E-06 

479 

1000.0 

331 

0.50E-05 

0.14E-06 

530 

1200.0 

395 

0.44E-05 

0.13E-06 

632 

1400.0 

459 

0.38E-05 

O.llE-06 

734 

1600.0 

523 

0.31E-05 

0.97E-07 

835 

1800.0 

587 

0.28E-05 

937 

2000.0 

651 

0.23E-05 

0.64E-07 

1038 

2400.0 

778 

0.29E-05 

0.15E-06 

1240 

2800.0 

906 

0.19E-05 

0.84E-07 

1442 

4000.0 

1288 

0.37E-05 

0.17E-06 

2047 
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Table  2:  Quadrature  performance  for  varying  precisions,  for  c  =  50 


e 

Maximum  Errors 

NpQi 

Roots 

Refined 

O.lOE-01 

19 

0.45E-01 

O.lOE-01 

30 

O.lOE-02 

20 

0.70E-02 

0.13E-02 

32 

O.lOE-03 

21 

0.91E-03 

0.14E-03 

33 

O.lOE-04 

22 

0.82E-04 

0.13E-04 

34 

O.lOE-05 

23 

0.54E-04 

O.llE-05 

36 

O.lOE-06 

24 

0.35E-05 

0.83E-07 

37 

O.lOE-07 

25 

0.33E-05 

0.57E-08 

38 

O.lOE-08 

26 

0.18E-06 

0.36E-09 

39 

O.lOE-09 

26 

0.18E-06 

0.36E-09 

40 

O.lOE-10 

27 

0.17E-06 

0.21E-10 

42 

O.lOE-11 

28 

0.79E-08 

O.llE-11 

43 

O.lOE-12 

29 

0.78E-08 

0.56E-13 

45 

O.lOE-13 

30 

0.31E-09 

0.27E-14 

55 

55 


Table  3:  Interpolation  performance  for  varying  band  limits,  for  e  =  10“’’ 


c 

Maximum  Errors 

NpoiX 

Roots 

Leg. 

5.0 

13 

0.12E-06 

0.12E-06 

17 

10.0 

18 

0.12E-06 

0.13E-06 

25 

15.0 

22 

0.24E-06 

0.25E-06 

31 

32 

20.0 

26 

0.26E-06 

0.28E-06 

37 

39 

25.0 

30 

0.22E-06 

0.23E-06 

43 

45 

30.0 

33 

0.67E-06 

0.73E-06 

49 

51 

35.0 

37 

0.42E-06 

0.46E-06 

55 

57 

40.0 

41 

0.25E-06 

0.27E-06 

61 

63 

45.0 

44 

0.54E-06 

0.60E-06 

67 

69 

50.0 

48 

0.29E-06 

0.33E-06 

73 

75 

100.0 

82 

0.39E-06 

0.46E-06 

128 

131 

150.0 

115 

0.52E-06 

0.64E-06 

182 

186 

200.0 

147 

0.12E-05 

0.15E-05 

235 

239 

250.0 

180 

0.83E-06 

O.llE-05 

287 

292 

300.0 

212 

0.13E-05 

0.17E-05 

340 

345 

350.0 

245 

0.75E-06 

O.lOE-05 

392 

398 

400.0 

277 

O.lOE-05 

0.14E-05 

443 

450 

450.0 

309 

0.13E-05 

0.18E-05 

495 

502 

500.0 

341 

0.16E-05 

0.22E-05 

547 

554 

1000.0 

662 

0.16E-05 

0.24E-05 

1058 

1068 

1500.0 

982 

0.15E-05 

0.25E-05 

1566 

1578 

2000.0 

1301 

0.20E-05 

0.35E-05 

2072 

2086 
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Table  4:  Interpolation  performance  for  varying  precisions,  for  c  =  25 


£ 

n 

Maximum  Errors 

■Npoi 

Roots 

Refined 

Cheb. 

Leg. 

O.lOE-01 

21 

0.38E-01 

0.43E-01 

31 

34 

O.lOE-02 

23 

0.37E-02 

0.41E-02 

34 

36 

O.lOE-03 

25 

0.29E-03 

0.31E-03 

37 

39 

O.lOE-04 

26 

0.74E-04 

0.81E-04 

39 

41 

O.lOE-05 

28 

0.44E-05 

0.47E-05 

41 

43 

O.lOE-06 

30 

0.22E-06 

0.23E-06 

43 

45 

O.lOE-07 

31 

0.46E-07 

0.49E-07 

45 

47 

O.lOE-08 

32 

0.95E-08 

O.lOE-07 

47 

49 

O.lOE-09 

34 

0.36E-09 

0.38E-09 

49 

51 

O.lOE-10 

35 

0.67E-10 

0.70E-10 

51 

52 

O.lOE-11 

37 

0.21E-11 

0.22E-11 

53 

54 

O.lOE-12 

38 

0.36E-12 

0.37E-12 

54 

56 

O.lOE-13 

39 

0.59E-13 

0.63E-13 

98 

61 

Table  5:  Quadrature  nodes  for  band-limited  functions,  with  c  =  50  and  £■  =  10  ^ 

This  table  contains  only  half  of  the  nodes  and  weights,  in  particular  those  for  which  the 
node  is  less  than  or  equal  to  zero;  reflecting  these  nodes  around  zero  yields  the  remaining 
nodes,  the  weight  for  the  node  at  —x  being  the  same  as  the  weight  for  the  node  at  x. 


Weight 


Node 

-.9904522459960804E+00 

-.9525601106643832E+00 

-.8927960861459153E+00 

-.8186117530609125E+00 

-.7350624131965875E+00 

-.6452878027260844E-f-00 

-.5512554698695428E+00 

-.4542505281525226E+00 

-.3551568458127944E+00 

-.2546173463813596E+00 

-.1531287781860989E-t-00 

-.5110121484050418E-01 


0.2413064234922188E-01 
0.5024347217095568E-01 
0.6801787677830858E-01 
0.7952155999100788E-01 
0.8706680708376023E-01 
0.9216240765763570E-01 
0.9569254015486106E-01 
0.9817257766311556E-01 
0.9990914516102242E-01 
0. 101088017264871 5E-f  00 
0.1018214308931439E4-00 
0.1021735189986602E-F00 
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Figure  2;  Maximum  norms  of  derivatives  of  prolate  spheroidal  wave  functions  for  c  =  200, 
and  of  normalized  Legendre  polynomials 
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Table  6:  Quadrature  nodes  for  band-limited  functions,  with  c  =  150  and  e  =  10 

This  table  contains  only  half  of  the  nodes  and  weights,  in  particular  those  for  which  the 
node  is  less  than  or  equal  to  zero;  reflecting  these  nodes  around  zero  yields  the  remaining 
nodes,  the  weight  for  the  node  at  —x  being  the  same  as  the  weight  for  the  node  at  x. 


Node 

-.9982883010959975E+00 

-.9911354691596528E+00 

-.9788315280982487E-I-00 

-.9621348937901911E+00 

-  .9418386698454396E-f  00 
-.9186509576802944E+00 
-.8931541850293142E+00 
-.8658083894041821E+00 

-  .8369709588254746E+00 
-.8069187108185302E-f00 
- .  7758670331396409E-t-00 
-.7439849501152674E+00 
- .  71 14064976175457E4-00 

-  .6782391686910609E-f 00 
-.6445701594098660E-h00 
-.6104710013384929E+00 

-  .5760010202980960E-h00 
-.5412099413257457E+00 
-.5061398697742787E4-00 

-  .4708268134473433E-1-00 
-.4353018643598344E+00 

-  .3995921259242572E+00 
-.3637214481257228E+00 
-.3277110167114320E-f-00 
-.2915798305819667E+00 

-  .2553450930388687E-f-00 
-.2190225363501577E-F00 
- .  1826266945721476E-F00 
- .  146171 1362450572E+00 
- .  1096686661347072E+00 
-.7313150339365902E-01 
-.3657144220122915E-01 
0 


Weight 

0.4374483371752129E-02 

0.9842619236149078E-02 

0.1463518300250369E-01 

0.1862396111287527E-01 

0.2184988739217138E-01 

0.2442858670932862E-01 

0.2648864579258096E-01 

0.2814375940413615E-01 

0.2948528624795690E-01 

0.3058356160435090E-01 

0.3149181066633766E-01 

0.3225015506203403E-01 

0.3288893713079314E-01 

0.3343126421620424E-01 

0.3389488931551181E-01 

0.3429358206877410E-01 

0.3463812513892117E-01 

0.3493704033879884E-01 

0.3519712095895683E-01 

0.3542382499917732E-01 

0.3562156808557525E-01 

0.3579394352776868E-01 

0.3594388900778062E-01 

0.3607381381247460E-01 

0.3618569660385742E-01 

0.3628116095737887E-01 

0.3636153393399723E-01 

0.3642789154364812E-01 

0.3648109393796617E-01 

0.3652181242257066E-01 

0.3655054982303338E-01 

0.3656765531685031E-01 

0.3657333451556860E-01 
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